Time Series Saham NVIDIA dengan Metode ARIMA & ARCH/GARCh

Angga Fathan Rofiqy

30 November, 2023

Kode di Hide dalam default, untuk menampilkan kode, klik Code .

Analisis sebelum ini : https://rpubs.com/ZenR_Prog/MPDW-Prak10

#                      -=( Install & Load Package Function )=-
install_load <- function (package1, ...)  {   

   # convert arguments to vector
   packages <- c(package1, ...)

   # start loop to determine if each package is installed
   for(package in packages){

       # if package is installed locally, load
       if(package %in% rownames(installed.packages()))
          do.call('library', list(package))

       # if package is not installed locally, download, then load
       else {
          install.packages(package)
          do.call("library", list(package))
       }
   } 
}

#Path Function
path <- function(){
  gsub  ( "\\\\",  "/",  readClipboard ()  )
}
#Copy path, Panggil function di console
#Copy r path, paste ke var yang diinginkan
#Export chart
export.chart <- "C:/Users/Fathan/Documents/Obsidian Vault/2. Kuliah/Smt 5/6. Metode Peramalan Deret Waktu/@Proj/STA1341-MPDW/Pertemuan 10-14/Chart2"

1 Pendahuluan

1.1 NVIDIA

Sumber : finance.yahoo.com

Terdapat beberapa faktor yang menjadi landasan dan membuat kami memilih saham NVIDA. Berikut adalah latar belakang yang melandasi pemilihan ini:

  1. Dominasi di Industri Chip AI
    NVIDIA Corporation tidak hanya merupakan pemimpin, tetapi juga menggambarkan dominasi dalam industri chip kecerdasan buatan (AI). Sejak pendiriannya pada tahun 1993, perusahaan ini telah memainkan peran kunci dalam menghadirkan solusi grafis dan komputasi untuk berbagai industri, termasuk gaming, visualisasi profesional, pusat data, dan otomotif.
  2. Inovasi Chip H200 sebagai Pemacu Utama
    Pengumuman pengembangan chip terbaru, H200, menambahkan dimensi baru pada potensi inovatif NVIDIA dalam mendukung perkembangan AI. Dengan peningkatan signifikan pada memori dan kecepatan, H200 memberikan daya ungkit yang kuat bagi aplikasi AI generatif dan model bahasa besar, menjadikannya subjek yang menarik untuk diteliti.
  3. Peran Penting dalam Layanan AI Generatif
    NVIDIA tidak hanya mendominasi pasar chip AI, tetapi juga memainkan peran penting dalam layanan AI generatif, termasuk layanan ChatGPT yang populer. Sebagai penyedia teknologi di balik berbagai aplikasi AI canggih, perusahaan ini memiliki dampak yang signifikan pada perkembangan dan kemajuan kecerdasan buatan.
  4. Dukungan dari Cloud Service Providers Terkemuka
    Keterlibatan NVIDIA dengan penyedia Cloud Service terkemuka seperti Amazon Web Services, Google Cloud, Microsoft Azure, dan Oracle Cloud Infrastructure menunjukkan pengakuan industri terhadap kualitas dan relevansi produk-produknya. Ini menambah nilai pada data saham NVIDIA sebagai subjek penelitian.
  5. Performa Saham yang Mengesankan
    Kinerja saham NVIDIA yang mengesankan, terutama dengan kenaikan lebih dari tiga kali lipat dalam nilai selama tahun ini, menarik perhatian sebagai indikator potensial untuk menganalisis dampak inovasi teknologinya terhadap nilai pasar saham.

Dengan kombinasi faktor-faktor ini, pemilihan data saham NVIDIA bukan hanya memungkinkan analisis yang mendalam tentang hubungan antara inovasi teknologi dan pergerakan harga saham, tetapi juga memberikan wawasan yang berharga tentang bagaimana NVIDIA terus memainkan peran krusial dalam perkembangan industri kecerdasan buatan.

1.2 Dataset

Dataset yang saya gunakan merupakan koleksi data harga saham historis periode 1 Januari 2018 hingga 10 November 2023 dari 6 saham yakni saham saingan NVDA seperti AMD, ARM, AVGO, TSM, dan INTC.

Dataset ini memilki data :

  1. Open: yakni Harga saham pada awal periode perdagangan tertentu. Ini adalah harga saham pertama pada hari perdagangan tersebut.
  2. High: Harga tertinggi yang saham capai selama periode perdagangan tersebut. Ini mencerminkan harga tertinggi yang pembeli bersedia bayar selama hari tersebut.
  3. Low: Harga terendah yang saham capai selama periode perdagangan tersebut. Ini mencerminkan harga terendah yang penjual bersedia terima selama hari tersebut.
  4. Close: Harga saham pada akhir periode perdagangan tertentu. Ini adalah harga saham terakhir pada hari perdagangan tersebut.
  5. Adj Close (Adjusted Close): Harga penutup yang telah disesuaikan untuk memperhitungkan perubahan seperti pembagian saham atau dividen. Ini adalah harga penutup yang paling relevan untuk analisis jangka panjang, karena mencerminkan harga saham yang sebenarnya setelah penyesuaian.
  6. Volume: Volume perdagangan saham selama periode tertentu. Ini mencerminkan jumlah saham yang diperdagangkan selama hari perdagangan tersebut.

Kami akan menggunakan peubah Adj Close (Adjusted Close), Karena sesuai dengan penjelasan diatas, peubah Adj Close adalah yang paling sesuai untuk dianalisis dibandingkan peubah lainnya.

1.3 Tujuan

Tujuan dari praktikum ini adalah untuk menganalisis pola perkiraan pergerakan harga tujuh saham teknologi terkemuka dengan harapan dapat memberikan rekomendasi kepada pembaca mengenai saham mana yang sebaiknya dipertimbangkan untuk dibeli atau diinvestasikan secara signifikan di antara tujuh perusahaan teknologi besar tersebut.

1.4 Data Preparation

1.4.1 Import Data

install_load('rio')
raw.data <- import("https://raw.githubusercontent.com/Zen-Rofiqy/STA1341-MPDW/main/Data/New/%40AAANTI%20Stock%20Prices.csv")

1.4.2 Data Checking

Cek Tipe data.

str(raw.data)
## 'data.frame':    7464 obs. of  8 variables:
##  $ Name     : chr  "AMD" "AMD" "AMD" "AMD" ...
##  $ Date     : IDate, format: "2018-01-02" "2018-01-03" ...
##  $ Open     : num  10.4 11.6 12.1 12.2 12 ...
##  $ High     : num  11 12.1 12.4 12.2 12.3 ...
##  $ Low      : num  10.3 11.4 12 11.7 11.8 ...
##  $ Close    : num  11 11.6 12.1 11.9 12.3 ...
##  $ Adj Close: num  11 11.6 12.1 11.9 12.3 ...
##  $ Volume   : int  44146300 154066700 109503000 63808900 63346000 62560900 52561200 38354900 47149300 42686600 ...

Semua data Karakter, harus diubah.

Cek Data kosong.

sum(is.na(raw.data))
## [1] 0

Tidak ada data kosong.

1.4.3 Penyesuaian Tipe Data

Semua tipe data masih berupa character. Harus diubah menjadi tipe data yang sesuai.

install_load('dplyr')
data <- raw.data %>%  
  mutate(
    Date = as.Date(raw.data[, 2], format = "%m/%d/%y"), #Mengubah menjadi Date 
    across(3:ncol(raw.data), as.numeric)                #Mengubah menjadi Numerik
  ) 
str(data)
## 'data.frame':    7464 obs. of  8 variables:
##  $ Name     : chr  "AMD" "AMD" "AMD" "AMD" ...
##  $ Date     : Date, format: "2018-01-02" "2018-01-03" ...
##  $ Open     : num  10.4 11.6 12.1 12.2 12 ...
##  $ High     : num  11 12.1 12.4 12.2 12.3 ...
##  $ Low      : num  10.3 11.4 12 11.7 11.8 ...
##  $ Close    : num  11 11.6 12.1 11.9 12.3 ...
##  $ Adj Close: num  11 11.6 12.1 11.9 12.3 ...
##  $ Volume   : num  4.41e+07 1.54e+08 1.10e+08 6.38e+07 6.33e+07 ...

1.4.4 Rechecking Data 

Cek kembali data kosong.

cat('Banyaknya Data Kosong', sum(is.na(data)))
## Banyaknya Data Kosong 0

1.4.5 Cek Periode Data

data2 <- data
install_load("lubridate")

dates <- as.Date(data2$Date)

# Buat rentang waktu mulai dari tanggal pertama hingga tanggal terakhir dalam data
full_date_range <- seq(min(dates), max(dates), by = "days")

# Bandingkan rentang waktu dengan tanggal yang ada dalam data
missing_dates <- setdiff(full_date_range, dates) 

# Jika 'missing_dates' kosong, maka semua tanggal sudah ada dalam data
if (length(missing_dates) == 0) {
  cat("Semua tanggal ada dalam data.\n")
} else {
  cat("Tanggal yang tidak ada dalam data sebanyak", length(missing_dates),
      "\nAtau sebanyak", length(missing_dates) * 6, 
      "Data Hilang dari ke-6 perusahaan yang ada")
}
## Tanggal yang tidak ada dalam data sebanyak 667 
## Atau sebanyak 4002 Data Hilang dari ke-6 perusahaan yang ada

Inputasi Data

install_load('purrr')
# Fungsi untuk mengisi data yang hilang
fill_missing_data <- function(name) {
  data_filtered <- data2 %>%
    filter(Name == name)
  
  full_date_range <- seq(min(data2$Date), max(data2$Date), by = "days")
  data_frame_template <- data.frame(Date = full_date_range)
  
  # Menambahkan kolom "Name" sesuai dengan perusahaan yang diproses
  data_frame_template$Name <- name
  
  data_filled <- merge(data_frame_template, data_filtered, 
                       by = c("Date", "Name"), all.x = TRUE)
  return(data_filled)
}

# Menggunakan purrr::map untuk memproses setiap nama perusahaan
filled_data_list <- map(unique(data2$Name), fill_missing_data)

# Gabungkan data-data yang telah diisi menjadi satu data frame
final_data <- data.frame()
final_data <- do.call(rbind, filled_data_list)

# Urutkan data berdasarkan "Name" terlebih dahulu, kemudian "Date"
final_data <- final_data %>%
  dplyr::select(1, 2, 7) %>%
  arrange(Name, Date)

#Input Data Hilang
install_load('imputeTS')

data <- na_interpolation(final_data$`Adj Close`, option = "linear")

install_load('ggplot2')
chart <- ggplot_na_imputations(final_data$`Adj Close`, data)
chart

#Export Chart
ggsave("00_Imputasi Data.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 23)

data <- final_data %>% dplyr::select(Name, Date) %>% 
  mutate(`Adj Close` = data)

Data sudah di imputasi.

Cek ukuran data

cat("Ukuran data awal adalah", nrow(data2), "(", nrow(data2)/6, "Jika persaham)",
    "\nBanyaknya peride data yang hilang", (nrow(data)-nrow(data2))/6,
    "(PerSaham.", (nrow(data)-nrow(data2)) , "Jika semua)",
    "\nUkuran data yang baru seharusnya =", 
    nrow(data2) + (nrow(data)-nrow(data2)) ,
    "\nIni sudah sesuai dengan Ukuran data input yakni =", nrow(data),
    "(", nrow(data)/6, "jika persaham)",
    "\nJadi ada sekitar", (nrow(data)-nrow(data2))/nrow(data)*100,"% Data hilang")
## Ukuran data awal adalah 7464 ( 1244 Jika persaham) 
## Banyaknya peride data yang hilang 906 (PerSaham. 5436 Jika semua) 
## Ukuran data yang baru seharusnya = 12900 
## Ini sudah sesuai dengan Ukuran data input yakni = 12900 ( 2150 jika persaham) 
## Jadi ada sekitar 42.13953 % Data hilang

Data sudah benar, siap untuk dianalisis.
Jika dalam rentang tahun 2021 hingga 2023 akan seperti :

coba <- data %>% 
  filter(Date >= as.Date("2021-01-01") & Date <= as.Date("2023-11-10")) 
coba2 <- data2 %>% 
  filter(Date >= as.Date("2021-01-01") & Date <= as.Date("2023-11-10")) 

cat("Ukuran data awal adalah", nrow(coba2), "(", nrow(coba2)/6, "Jika persaham)",
    "\nBanyaknya peride data yang hilang", (nrow(coba)-nrow(coba2))/6,
    "(PerSaham.", (nrow(coba)-nrow(coba2)) , "Jika semua)",
    "\nUkuran data yang baru seharusnya =", 
    nrow(coba2) + (nrow(coba)-nrow(coba2)) ,
    "\nIni sudah sesuai dengan Ukuran data input yakni =", nrow(coba),
    "(", nrow(coba)/6, "jika persaham)",
    "\nJadi ada sekitar", (nrow(coba)-nrow(coba2))/nrow(coba)*100,"% Data hilang")
## Ukuran data awal adalah 3642 ( 607 Jika persaham) 
## Banyaknya peride data yang hilang 437 (PerSaham. 2622 Jika semua) 
## Ukuran data yang baru seharusnya = 6264 
## Ini sudah sesuai dengan Ukuran data input yakni = 6264 ( 1044 jika persaham) 
## Jadi ada sekitar 41.85824 % Data hilang

1.4.6 Data Cleaned

install_load('DT')
datatable(data, filter = 'top', 
          options = list(pageLength = 5))

2 Eksplorasi Data

Referensi : Warna1, Warna2

col.amd <- c("#D02A49", "#F3AEAA"); col.arm <- c("#4EC2C1", "#B8E0DC") 
col.avgo <- c("#DC4F26", "#FCCFBC"); col.nvda <- c("#36B450", "#ACD694") 
col.tsm <- c("#EDB64E", "#F6DEB3"); col.intc <- c("#4B75BA", "#9BC6EA") 

cols <- c("AMD"=col.amd[1],   "test_amd"=col.amd[2],
          "ARM"=col.arm[1],   "test_arm"=col.arm[2],
          "AVGO"=col.avgo[1], "test_avgo"=col.avgo[2],
          "NVDA"=col.nvda[1], "test_mvda"=col.nvda[2],
          "TSM"=col.tsm[1],   "test_tsm"=col.tsm[2],
          "INTC"=col.intc[1], "test_intc"=col.intc[2])

2.1 Plot Time Series

install_load('ggplot2','extrafont')
# font_import(); loadfonts() #Run ini sekali aja
theme.ts <- list(
  theme(legend.position = "none",
        axis.text.x = element_text(hjust = 1, 
                                   margin = margin(b = 10, t=20)),
        axis.text.y = element_text(vjust = 0.5, face = "bold", 
                                   margin = margin(l = 20, r = 20)),
        plot.title = element_text(hjust = 0.5, face = "bold"),
        text = element_text(size = 30),
        plot.subtitle = element_text(hjust = 0.5),
        panel.background = element_rect(fill = 'transparent'),
        plot.background = element_rect(fill='transparent', color=NA),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(linewidth = 1, colour = "black"))
        )
theme.ts1 <- list(
  theme(legend.position = "none",
        axis.text.x = element_text(hjust = 1, 
                                   margin = margin(b = 10, t=20)),
        axis.text.y = element_text(vjust = 0.5, face = "bold", 
                                   margin = margin(l = 50, r = 20)),
        plot.title = element_text(hjust = 0.5, face = "bold"),
        text = element_text(size = 30),
        plot.subtitle = element_text(hjust = 0.5),
        panel.background = element_rect(fill = 'transparent'),
        plot.background = element_rect(fill='transparent', color=NA),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(linewidth = 1, colour = "black"))
        )

2.1.1 AAANTI

Melihat keseluruhan Time Series data saham.

install_load('viridis','ggrepel')
#Plot
chart <-
ggplot(data, aes(x=Date, y=`Adj Close`, color=Name, alpha=Name)) + #Data
  geom_line(aes(color=Name), linewidth=1.5) + #Timeseries
  #Color
  scale_color_manual(values = c(AMD=col.amd[1],   ARM=col.arm[1],
                                AVGO=col.avgo[1], NVDA=col.nvda[1],
                                TSM=col.tsm[1],   INTC=col.intc[1] )) +
  scale_alpha_manual(values = c("AMD" = .25, "ARM" = .25, "AVGO" = .25, 
                                "NVDA" = 1, "TSM" = .25, "INTC" = .25)) +
  theme.ts + #THeme
  labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
       title = "Time Series NVDA",
       subtitle = "Seperti apa sih pola deret waktu saham NVDA dan saingannya?\n") +
  # Label / legend
  geom_text_repel(
    data=data[data$Date == max(data$Date),], #Posisi di ujung data
    aes(color = Name, label = Name), #Warna garis & label saham
    size = 8, #Ukuran text
    nudge_x = 80, #Posisi Text (kanan 50)
    hjust = 0, #Ujung
    segment.size = 1,               #Ukuran garis
    segment.alpha = .75,             #transparasi garis
    segment.linetype = "dotted",    #Time garis
    box.padding = .4, #Biar label saham nggak dempetan
    segment.curvature = -0.1, #biar garis mulus
    segment.ncp = 8, 
    segment.angle = 60 
  ) +
  #Axis
    coord_cartesian(clip = "off"
  ) +
    scale_x_date( #Sumbu x
    date_breaks = "1 year",  # Menampilkan label setiap tahun
    date_labels = "%Y",  # Format label tahun
    limits = c(as.Date(min(data$Date)), 
               as.Date(max(data$Date)) + 120)
    #Tampilin lebih dari 20023-07-28 agar label saham bisa masuk
  ) +
    scale_y_continuous( #Sumbu y
    labels = scales::dollar_format(prefix = "$") #tambahin dolar
  ) +
    annotate( #Buat nandain batas data
    "text", x = as.Date(max(data$Date)), y = -35, 
    label = max(data$Date), size=6
  ) +
  geom_vline( #Buat garis batas data
    xintercept = as.numeric(as.Date( max(data$Date) )), 
             linetype = "dotted", color = "red")
chart

#Export Chart
ggsave("01_Time Series AAANTI.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 27)

Data saham berakhir pada tanggal 21 November 2023 dengan harga saham NVDA merupakan saham tertinggi kedua. Jika dilihat dari tahun 2019-2022, semua saham cenderung memiliki pola trend naik. Lalu dari 2022-2023 polanya cenderung trend turun. Dan 2023 keatas, polanya cenderung naik.

data2 <- data %>% 
  filter(Date >= as.Date("2021-01-01") & Date <= as.Date("2023-11-10")) 
install_load('viridis','ggrepel')
#Plot
chart <-
ggplot(data2, aes(x=Date, y=`Adj Close`, color=Name, alpha=Name)) + #Data
  geom_line(aes(color=Name), linewidth=1.5) + #Timeseries
  #Color
  scale_color_manual(values = c(AMD=col.amd[1],   ARM=col.arm[1],
                                AVGO=col.avgo[1], NVDA=col.nvda[1],
                                TSM=col.tsm[1],   INTC=col.intc[1] )) +
  scale_alpha_manual(values = c("AMD" = .35, "ARM" = .35, "AVGO" = .35, 
                                "NVDA" = 1, "TSM" = .35, "INTC" = .35)) +
  theme.ts + #THeme
  labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
       title = "Time Series NVDA 2021-2023",
       subtitle = "Seperti apa sih pola deret waktu saham NVDA dan saingannya?\n") +
  # Label / legend
  geom_text_repel(
    data=data2[data2$Date == max(data2$Date),], #Posisi di ujung data
    aes(color = Name, label = Name), #Warna garis & label saham
    size = 8, #Ukuran text
    nudge_x = 20, #Posisi Text (kanan 50)
    hjust = 0, #Ujung
    segment.size = 1,               #Ukuran garis
    segment.alpha = .75,             #transparasi garis
    segment.linetype = "dotted",    #Time garis
    box.padding = .4, #Biar label saham nggak dempetan
    segment.curvature = -0.1, #biar garis mulus
    segment.ncp = 8, 
    segment.angle = 60 
  ) +
  #Axis
    coord_cartesian(clip = "off"
  ) +
    scale_x_date( #Sumbu x
    date_breaks = "1 year",  # Menampilkan label setiap tahun
    date_labels = "%Y",  # Format label tahun
    limits = c(as.Date(min(data2$Date)), 
               as.Date(max(data2$Date)) + 60)
    #Tampilin lebih dari 20023-07-28 agar label saham bisa masuk
  ) +
    scale_y_continuous( #Sumbu y
    labels = scales::dollar_format(prefix = "$") #tambahin dolar
  ) +
    annotate( #Buat nandain batas data
    "text", x = as.Date(max(data2$Date)), y = -35, 
    label = max(data2$Date), size=6
  ) +
  geom_vline( #Buat garis batas data
    xintercept = as.numeric(as.Date( max(data2$Date) )), 
             linetype = "dotted", color = "red")
chart

#Export Chart
ggsave("01_Time Series AAANTI_2022-2023.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 27)

2.1.2 NVDA

nvda <- data2 %>%
  filter(Name == "NVDA")  # Filter data saham Amazon tahun 2022 ke atas

rownames(nvda) <- NULL
str(nvda)
## 'data.frame':    1044 obs. of  3 variables:
##  $ Name     : chr  "NVDA" "NVDA" "NVDA" "NVDA" ...
##  $ Date     : Date, format: "2021-01-01" "2021-01-02" ...
##  $ Adj Close: num  130 131 131 131 134 ...
datatable(nvda, filter = 'top', 
          options = list(pageLength = 5))

Mengubah Ajd Close Menjadi Time series.

nvda.ts <- ts(nvda[,3])

Ringkasan Data Ajd CLose.

summary(nvda.ts)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   112.2   158.1   206.4   237.1   278.5   493.5
install_load('ggtext')
min_value <- min(nvda$`Adj Close`)
min_date <- nvda$Date[which.min(nvda$`Adj Close`)]
min_perc <- (which.min(nvda$`Adj Close`) / nrow(nvda)) * 100

max_value <- max(nvda$`Adj Close`)
max_date <- nvda$Date[which.max(nvda$`Adj Close`)]
max_perc <- (which.max(nvda$`Adj Close`) / nrow(nvda)) * 100

chart <-
ggplot(nvda, aes(x=Date, y=`Adj Close`)) + 
  geom_line(aes(color=Name), linewidth=2) +
  scale_color_manual(values = col.nvda[1]) +
  labs(x = "\nPeriode (Tahun)", y='Saham Harga penutup',
       title = "Time Series Saham NVIDIA",
       subtitle = "Seperti apa sih pola deret waktu saham NVIDIA?\n") +
  theme(legend.position = "none") +
  theme.ts1 + 
  
#Ekstras

  #Titik terendah
  geom_segment(aes(x = min_date, 
                   xend = min_date, 
                   y = min(nvda$`Adj Close`)* 1.01, 
                   yend = max(nvda$`Adj Close`)*40/100), 
               arrow = arrow(type = "closed", length = unit(0.1, "inches")), 
               lineend = "round", color = "#D02A49", size=1.5) +
  geom_richtext(
    data = data.frame(x = min_date, 
                      y = max(nvda$`Adj Close`)*40/100, 
                      label = paste0("Titik Terendah") ),
    aes(x, y, label = label), size = 7, color = "white", 
    fill = "#D02A49", box.color = "white", parse = TRUE
  ) +
  geom_text(aes(x = min_date-1*35, y = max(`Adj Close`)*40/100, label = 
                  paste0(min_date)), 
            vjust = -1.5, hjust = 0, size = 7, color = "grey30") +

  #Titik tertinggi
  geom_segment(aes(x = max_date, 
                   xend = max_date, 
                   y = max(nvda$`Adj Close`)* .99, 
                   yend = max(nvda$`Adj Close`)*70/100), 
               arrow = arrow(type = "closed", length = unit(0.1, "inches")), 
               lineend = "round", color = "#4B75BA", size=1.5) +
  geom_richtext(
    data = data.frame(x = max_date, 
                      y = max(nvda$`Adj Close`)*70/100, 
                      label = paste0("Titik Tertinggi") ),
    aes(x, y, label = label), size = 7, color = "white", 
    fill = "#4B75BA", box.color = "white", parse = TRUE
  ) + 
  geom_text(aes(x = max_date-1*35, y = max(`Adj Close`)*60/100, label = 
                  paste0(max_date)), 
            vjust = -1.5, hjust = 0, size = 7, color = "grey30") +
  annotate( #Buat nandain batas data
    "text", x = as.Date(max(data2$Date)) -32, y = 50, 
    label = max(data2$Date), size=6
  ) +
  geom_vline( #Buat garis batas data
    xintercept = as.numeric(as.Date( max(data2$Date) )), 
             linetype = "dotted", color = "red", linewidth=1.5) 

chart

#Export Chart
ggsave("02_Time Series NVIDIA.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 27)

Berdasarkan grafik deret waktu yang disajikan, terlihat bahwa dari tahun 2021-2022 terjadi kenaikan harga saham/tren naik. Sedangkan 2022-2023 terjadi penurunan harga saham/tren turun. Dan pada tahun 2023 akhir tahun terjadi kenaikan harga saham yang drastis/tren naik.

2.1.3 Seasonal

install_load("forecast")
## Warning: package 'forecast' was built under R version 4.2.3
ggseasonplot(nvda.ts)
## Error in ggseasonplot(nvda.ts): Data are not seasonal

Terlihat bahwa data tidak memiliki pola musiman.

2.2 Data Train vs Test

#membagi 80% data latih (training) dan 20% data uji (testing)
train_nvda <- nvda[1: round(nrow(nvda) *78/100),]
test_nvda <- nvda[round(nrow(nvda) *78/100  +1): nrow(nvda),]
train_nvda.ts <- ts(train_nvda[,3])
test_nvda.ts <- ts(test_nvda[,3])
chart <-
ggplot() + 
  #Label Data Asli 
  annotate( "rect", alpha=.1, fill="#4EC2C1",
            xmin=as.Date(min(data2$Date)), 
            xmax=as.Date(data2$Date[nrow(train_nvda)]),
            ymin=max(nvda$`Adj Close`) * 1.25, ymax=Inf ) + 
  
  annotate( "text", color="#4EC2C1",
            x = as.Date(data2$Date[ceiling(nrow(train_nvda)/2)]), 
            y = max(nvda$`Adj Close`) * 1.3, 
    label = "Data Latih", size=10) + 
  annotate( "text", color="#4EC2C1",
            x = as.Date(data2$Date[ceiling(nrow(train_nvda)/2)]), 
            y = max(nvda$`Adj Close`) * 1.2, 
    label = paste0(nrow(train_nvda)," Hari (", 
                   round(nrow(train_nvda)/nrow(nvda)*100, 1), 
                   "%)"), size=10) + 
  
  #Label Data Ramal
  annotate( "rect", alpha=.1, fill="violetred",
            xmin=as.Date(data2$Date[nrow(train_nvda)]), 
            xmax=as.Date(max(data2$Date)),
            ymin=max(nvda$`Adj Close`) * 1.25, ymax=Inf ) + 
  
  annotate( "text", color="violetred",
            x = as.Date(data2$Date[ceiling( .89*length(data2$Date)/6)]) , 
            y = max(nvda$`Adj Close`) * 1.3, 
    label = "Data Uji", size=10) +
  annotate( "text", color="violetred",
            x = as.Date(data2$Date[ceiling( .89*length(data2$Date)/6)]) , 
            y = max(nvda$`Adj Close`) * 1.2, 
    label = paste0(nrow(test_nvda)," Hari (", 
                   round(nrow(test_nvda)/nrow(nvda)*100, 1), 
                   "%)"), size=10) +
  
  geom_vline( #Buat garis batas data
    xintercept = as.Date(data2$Date[nrow(train_nvda)]) , 
             linetype = "dotted", color = "black", linewidth = 2) +

#NVDA
  geom_line(data = train_nvda, linewidth=2,
            aes(x = Date, y = `Adj Close`, col = "NVDA")) +
  geom_line(data = test_nvda, linewidth=2,
            aes(x = Date, y = `Adj Close`, col = "test_nvda")) +
  geom_point(data = tail(train_nvda, 1), alpha = .5, 
             aes(x = Date, y = `Adj Close`), stroke=2,
             size = 15, shape = 21, color = "black", fill="violetred") +
  scale_colour_manual(values = c("NVDA"=col.nvda[1], "test_nvda"=col.nvda[2])) +
    theme.ts + #THeme
  labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
       title = "Time Series Saham NVDA",
       subtitle = "Pembagian Data Latih dan data Uji\n") +
  #Axis
    coord_cartesian(clip = "off"
  ) +
    scale_y_continuous( #Sumbu y
    labels = scales::dollar_format(prefix = "$") #tambahin dolar
  ) +

#Bagi Data Train & Test
  geom_segment(aes(x = as.Date(max(train_nvda$Date)), 
                   xend = as.Date(max(train_nvda$Date)), 
                   y = max(nvda$`Adj Close`)* .58, 
                   yend = max(nvda$`Adj Close`)*.7), 
               arrow = arrow(type = "closed", length = unit(0.1, "inches")), 
               lineend = "round", color = "#4B75BA", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(train_nvda$Date)), 
                      y = max(nvda$`Adj Close`)*.7, 
                      label = max(train_nvda$Date) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "#4B75BA", box.color = "white", parse = TRUE
  ) +

  geom_segment(aes(x = as.Date(max(data2$Date)), 
                   xend = as.Date(max(data2$Date)), 
                   y = max(nvda$`Adj Close`)* .95, 
                   yend = max(nvda$`Adj Close`)*1.05), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "#D02A49", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(data2$Date)), 
                      y = max(nvda$`Adj Close`)*1.05, 
                      label = max(data2$Date) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "#D02A49", box.color = "white", parse = TRUE
  )
  
chart

#Export Chart
ggsave("02_TS_NVDA_train-test.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 27)

Berdasarkan plot data deret waktu pada data latih (\(78\%\) dari data asli), terlihat bahwa data menunjukkan tren naik dan tren turun. Ini mengisyaratkan bahwa data latih tidak memenuhi kriteria stasioneritas dalam rataan maupun ragam. Di sisi lain, dalam plot data uji (\(22\%\) dari data asli), terlihat adanya tren yang melonjak naik dan kurangnya nilai tengah yang stabil. Ini juga menunjukkan bahwa data uji tidak stasioner dalam rataan.

3 Stasioneritas

3.1 Uji Stasioneritas

3.1.1 Plot ACF

install_load('tsibble','tseries')
acf(train_nvda.ts, main="", col=col.nvda[1], lwd=2)
mtext("NVDA", side=3, line=1, cex=2, font=2)

Berdasarkan plot ACF, terlihat bahwa plot ACF data train menurun secara perlahan (tails of slowly). Hal ini juga menjadi indikasi bahwa data tidak stasioner dalam rataan dan tidak membentuk gelombang sinus.

3.1.2 Uji ADF

tseries::adf.test(train_nvda.ts)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_nvda.ts
## Dickey-Fuller = -1.6692, Lag order = 9, p-value = 0.7183
## alternative hypothesis: stationary
if (tseries::adf.test(train_nvda.ts)[["p.value"]] < 0.05) {
  cat("Karena p-value < 0.05, Maka Rataan Stasioner")
} else {
  cat("Karena p-value > 0.05, Maka Rataan Tidak Stasioner")
}
## Karena p-value > 0.05, Maka Rataan Tidak Stasioner

\(H_0\) : Data tidak stasioner dalam rataan

\(H_1\) : Data stasioner dalam rataan

Berdasarkan uji ADF tersebut, p-value lebih besar dari taraf nyata \(5\%\) sehingga tak tolak \(H_0\) dan menandakan bahwa data tidak stasioner dalam rataan. Hal ini sesuai dengan hasil eksplorasi menggunakan plot time series dan plot ACF, sehingga ketidakstasioneran model kedepannya harus ditangani.

3.1.3 Plot Box-Cox

install_load('MASS')
index <- seq(1:nrow(train_nvda)) #Beda titik potong
bc =boxcox(train_nvda.ts~index, lambda=seq(-2, 4, by=.01))
title("NVDA", cex.main = 2)

lambda_nvda <- bc$x[which.max(bc$y)] #Nilai Rounded Lambda
sk_nvda <- bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)] #SK

Gambar di atas menunjukkan bahwa dari rentang nilai \(\lambda\) -2 hingga 4, pada selang yang dibuatnya tidak memuat nilai \(\lambda=1\), sehingga ragam tidak stasioner.

3.1.3.1 Lambda dan Selang Kepercayaan

lsk <- data.frame(
  Lambda = c(lambda_nvda ),
  Batas.Bawah = c( min(sk_nvda)  ),
  Batas.Atas = c(max(sk_nvda) )
) %>%
  mutate(across(where(is.numeric), ~round(., 2)))

lsk$Keterangan <- apply(lsk, 1, function(row) {
  ifelse(row["Batas.Bawah"] <= 1 && row["Batas.Atas"] >= 1,
         "Ragam Stasioner",  "Ragam Tidak Stasioner")
})

rownames(lsk) <- c("NVDA")

datatable(lsk, filter = "none")

Hal ini dibuktikan dengan tabel nilai rounded value (\(\lambda\)) dan selang kepercayaan diatas yang menyatakan bahwa memuat nilai satu di dalam selangnya. Ini mengindikasikan bahwa NVDA memerlukan transformasi data karena tidak memiliki sifat stasioner dalam ragam.

3.2 Penanganan Ketidakstasioneran Data

3.2.1 Dalam Ragam: Transformasi

Mengutip dari laman r-coder.com, formula transformasi data disesuaikan dengan nilai rounded value (\(\lambda\)) optimumnya.

Namun dari laman tersebut dan dari laman robjhyndman.com mengatakan bahwa rounded value (\(\lambda\)) = 1 itu cenderung tidak ada perubahan substansial atau signifikan yang terjadi pada data. Sehingga jika transformasi diatas tidak menyelesaikan masalah stasioneritas ragam maka alternatif transfromasi yang dapat dilakukan adalah dengan \((x^\lambda -1 )/\lambda\) dengan transformasi balik \((\lambda x+1)^{1/\lambda}\).
\(\mathbf{\lambda}\) Transformasi Transformasi Balik
-2 \(1/x^2\) \(1/\sqrt{x}\)
-1 \(1/x\) \(1/x\)
-0.5 \(1/\sqrt{x}\) \(1/x^2\)
0 \(\log{(x)}\) \(e^x\)
0.5 \(\sqrt{x}\) \(x^2\)
1 \(x\) \(x\)
2 \(x^2\) \(\sqrt{x}\)

Sehingga transformasi yang digunakan pada data saham NVDIA adalah sebagai berikut.

lamb <- c(-2, -1, -.5, 0, .5, 1, 2)
# Fungsi untuk mengonversi lambda ke transformasi
trans <- function(lambda) {
  lambda <- lamb[sapply(lambda, function(y) which.min(abs(y - lamb))) ]
  case_when(
    lambda == -2   ~ "1/x^2",
    lambda == -1   ~ "1/x",
    lambda == -0.5 ~ "1/sqrt(x)",
    lambda == 0    ~ "log(x)",
    lambda == 1    ~ "x",
    lambda == 0.5  ~ "sqrt(x)",
    lambda == 2    ~ "x^2",
    TRUE           ~ NA_character_
  )
}

# Fungsi untuk mengonversi lambda ke transformasi balik
trans_balik <- function(lambda) {
  lambda <- lamb[sapply(lambda, function(y) which.min(abs(y - lamb))) ]
  case_when(
    lambda == -2   ~ "1/sqrt(x)",
    lambda == -1   ~ "1/x",
    lambda == -0.5 ~ "1/x^2",
    lambda == 0    ~ "exp(x)",
    lambda == 1    ~ "x",
    lambda == 0.5  ~ "x^2",
    lambda == 2    ~ "sqrt(x)",
    TRUE           ~ NA_character_
  )
}

# Menambahkan kolom Transformasi dan Trans Balik
tran <- lsk %>%
  mutate(
    Transformasi = trans(Lambda),
    `Transformasi Balik` = trans_balik(Lambda)
  ) %>% dplyr::select(Lambda, Transformasi, `Transformasi Balik`)

datatable(tran, filter = 'none')

3.2.1.1 Transformasi

#Transformasi
train_nvda.ts.new <- 1/sqrt(train_nvda.ts)

#Plot BoxCox
index <- seq(1:nrow(train_nvda))
bc <-boxcox(train_nvda.ts.new ~index, lambda=seq(-2, 4, by=.01))
title("NVDA", cex.main=2)

lambda_nvda.new <- bc$x[which.max(bc$y)] #Nilai Rounded Lambda
sk_nvda.new <- bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)] #SK
lsk.new <- data.frame(
  Lambda = c(lambda_nvda.new ),
  Batas.Bawah = c( min(sk_nvda.new)  ),
  Batas.Atas = c( max(sk_nvda.new)  )
) %>%
  mutate(across(where(is.numeric), ~round(., 2)))

lsk.new$Keterangan <- apply(lsk.new, 1, function(row) {
  ifelse(row["Batas.Bawah"] <= 1 && row["Batas.Atas"] >= 1,
         "Ragam Stasioner",  "Ragam Tidak Stasioner")
})

rownames(lsk.new) <- c("NVDA")

datatable(lsk.new, filter = 'none')

Setelah melalui proses transformasi, gambar dan tabel di atas mengindikasikan bahwa setiap rangkaian data sekarang mencakup nilai satu dalam selangnya, yang menunjukkan bahwa data tersebut sekarang sudah stasioner dalam ragam.

3.2.1.2 Plot Transformasi

ragam <- data.frame(
  Date = train_nvda %>% dplyr::select(Date)  ,
  `Asli` = c(as.vector(train_nvda.ts) ),
  `Trans` = c(as.vector(train_nvda.ts.new) ) *3000,
  Name = rep(c("NVDA"), 
              each = nrow(train_nvda) )
)
colnames(ragam) <- c("Date", "Asli", "Trans", "Name")
chart <- 
ggplot() +
# Akhir data
  geom_vline( #Buat garis batas data
    xintercept = as.numeric(as.Date( max(train_nvda$Date) )), 
             linetype = "dotted", color = "red") +

#Time Series
  #Data Asli
  geom_line(data = ragam, linewidth=2,
            aes(x = Date, y = `Asli`, col = "Asli")) +
  #Data Transformasi
  geom_line(data = ragam, linewidth=2,
            aes(x = Date, y = `Trans`, col = "Trans")) +

  scale_colour_manual(values = c("Asli"=col.nvda[2], "Trans"=col.nvda[1])) +
  theme.ts + #THeme
  labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
       title = "Plot Transformasi NVDA",
       subtitle = "Perbandingan plot Transformasi dengan Plot Asli\n") +
  #Axis
    coord_cartesian(clip = "off"
  ) +
    scale_x_date( #Sumbu x
    date_breaks = "1 year",  # Menampilkan label setiap tahun
    date_labels = "%Y",  # Format label tahun
    limits = c(as.Date(min(ragam$Date)), 
               as.Date(max(ragam$Date)) + 80)
  ) +
    scale_y_continuous( #Sumbu y
    labels = scales::dollar_format(prefix = "$") #tambahin dolar
  ) +
  ylim(min(ragam$Asli), max(ragam$Asli)) +

#Data asli
  geom_segment(aes(x = as.Date(max(train_nvda$Date)), 
                   xend = as.Date(max(train_nvda$Date))+60, 
                   y = train_nvda$`Adj Close`[nrow(train_nvda)], 
                   yend = train_nvda$`Adj Close`[nrow(train_nvda)]), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "gray60", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(train_nvda$Date))+60, 
                      y = train_nvda$`Adj Close`[nrow(train_nvda)], 
                      label = "Data Asli" ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "gray60", box.color = "white", parse = TRUE
  ) +

# Data Transformasi
  geom_segment(aes(x = as.Date(max(train_nvda$Date)), 
                   xend = as.Date(max(train_nvda$Date)) + 60, 
                   y = ragam$Trans[nrow(ragam)], 
                   yend = ragam$Trans[nrow(ragam)] ), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "#D02A49", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(train_nvda$Date)) +60, 
                      y = ragam$Trans[nrow(ragam)], 
                      label = "Transformasi" ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "#D02A49", box.color = "white", parse = TRUE
  ) +
  
#Tanggal akhir
  geom_segment(aes(x = as.Date(max(train_nvda$Date)), 
                   xend = as.Date(max(train_nvda$Date)), 
                   y = -Inf, 
                   yend = max(train_nvda$`Adj Close`)*.15), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "gray60", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(train_nvda$Date)), 
                      y = max(train_nvda$`Adj Close`)*.15, 
                      label = max(train_nvda$Date) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "gray60", box.color = "white", parse = TRUE
  )

chart

#Export Chart
ggsave("03_Stas_4_Transform.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 27)

3.2.2 Dalam Rataan: Differencing

Differencing (diferensiasi) adalah suatu proses statistik yang digunakan untuk mengatasi ketidakstasioneran dalam data rata-rata atau tren data. Tujuannya adalah untuk membuat data menjadi lebih stasioner dengan mengurangi atau menghilangkan tren atau pola waktu yang mungkin ada dalam data tersebut.

#Diff
train_nvda.diff <- diff(train_nvda.ts.new, differences = 1) 

#Plot
plot.ts(train_nvda.diff, lty=1, xlab="Periode (Hari)", 
        col = col.nvda[1], lwd = 2,
        main="Plot Differencing", cex.main = 2)

Berdasarkan plot data deret waktu yang sudah di Differencing, terlihat bahwa seluruh data sudah stasioner dalam rataan ditandai dengan data bergerak pada nilai tengah tertentu (tidak terdapat trend ataupun musiman pada data).

3.3 Uji Ulang

Akan dilakukan uji ulang untuk mengecek kestasioneran data dalam rataan.

3.3.1 Plot ACF

stats::acf(train_nvda.diff, main="", col=col.nvda[1], lwd=2, ylim=c(-.1,.1))
mtext("NVDA", side=3, line=1, cex=2, font=2)

Berdasarkan ketujuh plot ACF tersebut, dapat dilihat bahwa lag nya sudah tidak tails off slowly (trennya telah berakhir dan korelasi antara pengamatan pada lag-lag yang lebih jauh telah berkurang secara signifikan), sehingga data dinyatakan stasioner dalam rataan. Ini menunjukkan bahwa masalah ketidakstasioneran dalam data telah berhasil diatasi.

3.3.2 Uji ADF

tseries::adf.test(train_nvda.diff)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_nvda.diff
## Dickey-Fuller = -8.3563, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary
if (tseries::adf.test(train_nvda.diff)[["p.value"]] < 0.05) {
  cat("Karena p-value < 0.05, Maka Rataan Stasioner")
} else {
  cat("Karena p-value > 0.05, Maka Rataan Tidak Stasioner")
}
## Karena p-value < 0.05, Maka Rataan Stasioner

\(H_0\) : Data tidak stasioner dalam rataan

\(H_1\) : Data stasioner dalam rataan

Berdasarkan uji ADF tersebut, p-value yang didapatkan yakni lebih kecil dari taraf nyata \(5\%\) sehingga tolak \(H_0\) atau data stasioner dalam rataan. Hal ini sesuai dengan hasil eksplorasi menggunakan plot time series dan plot ACF, sehingga dalam hal ini ketidakstasioneran data sudah berhasil ditangani dan dapat dilanjutkan ke pemodelan.

4 Model

4.1 Model Tentatif

Identifikasi model tentatif dapat dilakukan dengan melihat pada plot ACF, PACF dan Tabel EACF. Mengutip dari postingan pada laman medium.com Panduan dasar untuk menginterpretasikan plot ACF dan PACF adalah sebagai berikut :

  1. Cari pola tail off pada ACF atau PACF.
  2. Jika tail off terjadi pada ACF → model ARcut off pada PACF akan memberikan nilai p untuk AR(p).
  3. Jika tail off terjadi pada PACF → model MAcut off pada ACF akan memberikan nilai q untuk MA(q).
  4. Jika terjadi tail off pada kedua ACF dan PACF → model ARMA.

4.1.1 Plot ACF & PACF

par(mfrow=c(1,2))
stats::acf(train_nvda.diff, main="", col=col.nvda[1], lwd=2, ylim=c(-.1,.1))
mtext("NVDA", side=3, line=1, cex=2, font=2)

stats::pacf(train_nvda.diff, main="", col=col.nvda[1], lwd=2, ylim=c(-.1,.1))
mtext("NVDA", side=3, line=1, cex=2, font=2)

Referensi : Significance level of ACF and PACF in R

acf.pacf <- function(dt){ #Fungsi Deteksi Cut off plot ACF & PACF
  N <- acf(dt, plot = F)[["n.used"]] + 1 #Banyaknya Data
  #Garis Signifikan
  sl <- c(-1,1)*(exp(2*1.96/sqrt(N-3))-1)/(exp(2*1.96/sqrt(N-3))+1)
  acf <- data.frame( lag = acf(dt, plot=F)[["lag"]], 
                     acf = acf(dt, plot=F)[["acf"]] ) 
  
    cut.off1 <- acf %>%
        filter(acf <= min(sl) | acf >= max(sl)) %>%
        dplyr::select(lag); colnames(cut.off1) <- "acf"
  
  pacf <- data.frame(lag = pacf(dt, plot=F)[["lag"]], 
                     acf = pacf(dt, plot=F)[["acf"]] ) 
  
    cut.off2 <- pacf %>%
        filter(acf <= min(sl) | acf >= max(sl)) %>%
        dplyr::select(lag); colnames(cut.off2) <- "pacf"
    
  cut.off <- merge(cut.off1, cut.off2, by = 0, all = TRUE)[-1]
  return(cut.off)
}

data_frames <- list(
  train_nvda.diff
)

name <- "NVDA"

acf.pacf.tab <- data.frame()
for(i in 1:length(name)){
  dex <- data_frames[[i]] %>% acf.pacf() %>% 
  rename(!!paste0(name[i], ".a") := acf, 
         !!paste0(name[i], ".p") := pacf)
  
    acf.pacf.tab <- merge(acf.pacf.tab, dex, by=0, all=TRUE)[-1]
}

datatable(acf.pacf.tab, filter = 'none')

Berdasarkan plot tersebut, terlihat bahwa plot ACF cut off pada lag 0. Sedangkan pada plot PACF cut off pada lag 14. Sehingga model tidak dapat di identifikasi dengan plot ACF.

4.1.2 Tabel EACF

Referensi : [University of South Carolina] Chapter 6: Model Specification for Time Series

Mengutip dari referensi diatas, Identifikasi model menggunakan Tabel EACF (Extended Autocorrelation Function) untuk proses ARMA(p, q) seharusnya secara teoritis memiliki pola segitiga nol, dengan nol di sudut kiri atas muncul pada baris ke-p dan kolom ke-q (dengan label baris dan kolom dimulai dari 0). Pola segitiga ini merupakan pola segitiga pertama dari setiap baris.

Sebagai contoh :

Dalam hal ini model tentatif yang terbentuk adalah ARIMA(0,1,1), ARIMA(1,1,1), ARIMA(2,1,2), ARIMA(3,1,3), dst.. Namun karena kemungkinannya sangat banyak, maka akan digunakan function agar menyingkat proses.

Pemilihan model ARIMA terbaik dilakukan ketika memiliki nilai AIC yang terkecil, Semua Parameter Signifikan, dan Tidak ada Parameter NA. Kami akan mencoba untuk menggunakan function pemodelam model ARIMA berdasarkan tabel EACF. Namun dengan melakukan 2 pemodelan perbaris AR(p) nya, agar sekalian overfitting.

install_load('TSA', 'forecast')

# Fungsi untuk menghitung model ARIMA dan menganalisis parameter
model_tentatif <- function(data, p_max, d, q_max, alpha=0.05, drift=FALSE) {
  best_model <- NULL
  best_aic <- Inf
  eacf_result <- eacf(data) #Untuk Ektraksi & identifikasi model
  models <- data.frame(Model = character(0), 
                       AIC = numeric(0), 
                       Signif = character(0), 
                       Keterangan = character(0))
  
  for (p in 0:p_max) {
    overfit <- 0
    for (q in 1:q_max) {
      #Pola Matriks segitiga bawah
      if (!is.na(eacf_result$symbol[p + 1, q ]) && 
          !is.na(eacf_result$symbol[p + 1, q + 1]) && 
          !is.na(eacf_result$symbol[p + 2, q + 1])) {
        if (eacf_result$symbol[p + 1, q ] == "o" && 
            eacf_result$symbol[p + 1, q + 1] == "o" && 
            eacf_result$symbol[p + 2, q + 1] == "o") {
      
          model <- Arima(data, order=c(p,d,q), method="ML", 
                         include.drift=drift)
          aic <- AIC(model)
          
          # Mendapatkan nilai coef dari model
          coeftest_result <- lmtest::coeftest(model)
          
          # jika lebih kecil dari alpha, maka signifikan
          significant_params <- 
            rownames(coeftest_result)[coeftest_result[, "Pr(>|z|)"] < alpha]  
          
          # jika lebih besar dari alpha, maka tidak signifikan
          non_significant_params <- 
            rownames(coeftest_result)[coeftest_result[, "Pr(>|z|)"] > alpha]  
          
          # Keterangan signifikansi
          if (length(significant_params) == 0) {
            keterangan <- "Semua parameter tidak signifikan"
          } else if (length(significant_params) == nrow(coeftest_result)) {
            keterangan <- "Semua parameter signifikan"
          } else {
            keterangan <- paste("Parameter yang tidak signifikan adalah", 
                                paste(non_significant_params, collapse = ", "))
          }
          
          models <- rbind(models, 
                    data.frame(Model = paste("ARIMA(", p, ",", d, ",", q, ")", 
                                             sep = ""), 
                               AIC = aic, 
                               Signif = paste(significant_params, 
                                              collapse = ", "), 
                               Keterangan = keterangan))
          
          #Identifikasi Best Model
            if (keterangan == "Semua parameter signifikan" && 
                !any(is.na(significant_params))) {
              if (aic < best_aic) {
                best_model <- model
                best_aic <- aic
            }
          }
          #jika model ditemukan pada baris p. Maka lanjut ke p+1
          if(overfit >= 1) break
          overfit = overfit + 1
        }
      }
    }
  }
  
  cat("\nModel ARIMA dengan AIC terkecil:\n")
  print(best_model)

  return(models)
}
model.tentaif_nvda <- 
  model_tentatif(train_nvda.diff, p_max = 6, d = 1, q_max = 12)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o o o o o  o  o  o 
## 1 x o o o o o o o o o o  o  o  o 
## 2 x x o o o o o o o o o  o  o  o 
## 3 x o x o o o o o o o o  o  o  o 
## 4 o x x x o o o o o o o  o  o  o 
## 5 x x x o x o o o o o o  o  o  o 
## 6 x x x x x x o o o o o  o  o  o 
## 7 x x o x x x o o o o o  o  o  o 
## 
## Model ARIMA dengan AIC terkecil:
## Series: data 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -0.9933
## s.e.   0.0078
## 
## sigma^2 = 9.868e-07:  log likelihood = 4460.68
## AIC=-8917.36   AICc=-8917.34   BIC=-8907.96
datatable(model.tentaif_nvda, filter = 'top', 
          options = list(pageLength = 7))

Berdasarkan pendugaan parameter di atas, nilai AIC terkecil \(-8917.35\) dimiliki oleh model ARIMA(0,1,1) dan seluruh parameternya signifikan sehingga model yang dipilih adalah model ARIMA(0,1,1).

model_nvda.da <- Arima(train_nvda.diff, order=c(0,1,1), method="ML")

4.2 Model ARIMA dengan Drift

Tren menentukan bagaimana perubahan deret waktu dari waktu ke waktu. Dalam model ARIMA, ini dapat dimodelkan menggunakan persamaan dasar berikut “a” disebut “intercept” dan “b” disebut “drift”. “Drift” adalah hanyalah kemiringan garis lurus. Selanjutnya akan dilakukan kombinasi model ARIMA dengan Intercept dan Drift.

4.2.1 1. Model Tanpa Intercept dan Tanpa Drift

Simpelnya ini merupakan model data train yang diberi differencing. Yakni memodelkan data train dengan ordo I dari ARIMA(p, d, q) sama dengan 1, sehingga kan membentuk model ARIMA(p, 1, q).

drift1_nvda <- 
  model_tentatif(train_nvda.ts, p_max = 6, d = 1, q_max = 12)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x x x x  x  x  x 
## 1 o o o o o o o o o o o  o  o  o 
## 2 x o o o o o o o o o o  o  o  o 
## 3 x x o o o o o o o o o  o  o  o 
## 4 x x x o o o o o o o o  o  o  o 
## 5 x x x o o o o o o o o  o  o  o 
## 6 x x o x x o o o o o o  o  o  o 
## 7 x x x x x x o o o o o  o  o  o 
## 
## Model ARIMA dengan AIC terkecil:
## Series: data 
## ARIMA(4,1,4) 
## 
## Coefficients:
##          ar1      ar2     ar3      ar4      ma1     ma2      ma3     ma4
##       0.4960  -0.7914  0.5091  -0.9828  -0.4804  0.7803  -0.5046  0.9699
## s.e.  0.0144   0.0120  0.0094   0.0172   0.0236  0.0193   0.0172  0.0240
## 
## sigma^2 = 29.16:  log likelihood = -2522.64
## AIC=5063.29   AICc=5063.51   BIC=5105.59
datatable(drift1_nvda, filter = 'top', 
          options = list(pageLength = 7))

Berdasarkan pendugaan parameter di atas, nilai AIC terkecil \(5063.53\) dimiliki oleh model ARIMA(4,1,4) dan seluruh parameternya signifikan sehingga model yang dipilih adalah model ARIMA(4,1,4).

model_nvda.drift1 <- Arima(train_nvda.ts, order=c(4,1,4), method="ML")

4.2.2 2. Model Dengan Intercept Tanpa Drift

Simpelnya ini merupakan model data differencing tanpa diberi differencing. Yakni memodelkan data differencing dengan ordo I dari ARIMA(p, d, q) sama dengan 0, sehingga kan membentuk model ARIMA(p, 0, q).

drift2_nvda <- 
  model_tentatif(train_nvda.diff, p_max = 6, d = 0, q_max = 12)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o o o o o  o  o  o 
## 1 x o o o o o o o o o o  o  o  o 
## 2 x x o o o o o o o o o  o  o  o 
## 3 x o x o o o o o o o o  o  o  o 
## 4 o x x x o o o o o o o  o  o  o 
## 5 x x x o x o o o o o o  o  o  o 
## 6 x x x x x x o o o o o  o  o  o 
## 7 x x o x x x o o o o o  o  o  o 
## 
## Model ARIMA dengan AIC terkecil:
## NULL
datatable(drift2_nvda, filter = 'top', 
          options = list(pageLength = 7))

Kategori Model terbaik tidak ada yang terpenuhi, yakni tidak ada nilai AIC terkecil dengan semua parameter nya signifikan dan tidak ada NA di dalam nya. Sehingga Model Dengan Intercept Tanpa Drift pada NVDA tidak ada.

model_nvda.drift2 <- NA

4.2.3 3. Model Tanpa Intercept Dengan Drift

Simpelnya ini merupakan model data train yang diberi differencing dan drift. Yakni memodelkan data train dengan ordo I dari ARIMA(p, d, q) sama dengan 1, sehingga kan membentuk model ARIMA(p, 1, q) + include.drift.

drift3_nvda <- 
  model_tentatif(train_nvda.ts, p_max = 6, d = 1, q_max = 12, drift=TRUE)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x x x x  x  x  x 
## 1 o o o o o o o o o o o  o  o  o 
## 2 x o o o o o o o o o o  o  o  o 
## 3 x x o o o o o o o o o  o  o  o 
## 4 x x x o o o o o o o o  o  o  o 
## 5 x x x o o o o o o o o  o  o  o 
## 6 x x o x x o o o o o o  o  o  o 
## 7 x x x x x x o o o o o  o  o  o 
## 
## Model ARIMA dengan AIC terkecil:
## NULL
datatable(drift3_nvda, filter = 'top', 
          options = list(pageLength = 7))

Kategori Model terbaik tidak ada yang terpenuhi, yakni tidak ada nilai AIC terkecil dengan semua parameter nya signifikan dan tidak ada NA di dalam nya. Sehingga Model Tanpa Intercept Dengan Drift pada NVDA tidak ada.

model_nvda.drift3 <- NA

4.2.4 4. Model dengan Intercept dan dengan Drift

Simpelnya ini merupakan model data differencing tanpa diberi differencing namun diberi drift. Yakni memodelkan data differencing dengan ordo I dari ARIMA(p, d, q) sama dengan 0, sehingga kan membentuk model ARIMA(p, 0, q) + include.drift.

drift4_nvda <- 
  model_tentatif(train_nvda.diff, p_max = 6, d = 1, q_max = 12, drift=TRUE)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o o o o o  o  o  o 
## 1 x o o o o o o o o o o  o  o  o 
## 2 x x o o o o o o o o o  o  o  o 
## 3 x o x o o o o o o o o  o  o  o 
## 4 o x x x o o o o o o o  o  o  o 
## 5 x x x o x o o o o o o  o  o  o 
## 6 x x x x x x o o o o o  o  o  o 
## 7 x x o x x x o o o o o  o  o  o 
## 
## Model ARIMA dengan AIC terkecil:
## NULL
datatable(drift4_nvda, filter = 'top', 
          options = list(pageLength = 7))

Kategori Model terbaik tidak ada yang terpenuhi, yakni tidak ada nilai AIC terkecil dengan semua parameter nya signifikan dan tidak ada NA di dalam nya. Sehingga Model dengan Intercept dan dengan Drift pada NVDA tidak ada.

model_nvda.drift4 <- NA

4.3 Model Terbaik

model_terbaik <- function(model, rowname, mods) {
  if (!all(is.na(model))) {
    model_df <- data.frame(
      Model = paste("ARIMA(", 
                     paste(as.character(model[["call"]][["order"]][-1]), 
                           collapse = ","), ")", sep = ""),
      AIC = model[["aicc"]],
      row.names = rowname
    )
    model_df$Keterangan <- 
      mods[which(mods$Model == 
                     model_df["Model"][rowname,]), "Keterangan"]
  return(model_df)
    
  } else {
    return(data.frame(
      Model = NA_character_,
      AIC = NA_real_,
      Keterangan = NA_character_,
      row.names = rowname
    ))
  }
}

best_model <- rbind(
  model_terbaik(model_nvda.da, 'NVDA', model.tentaif_nvda),
  model_terbaik(model_nvda.drift1, 'NVDA', drift1_nvda),
  model_terbaik(model_nvda.drift2, 'NVDA', drift2_nvda),
  model_terbaik(model_nvda.drift3, 'NVDA', drift3_nvda),
  model_terbaik(model_nvda.drift4, 'NVDA', drift4_nvda)
) %>%
  mutate(across(where(is.numeric), ~round(., 2)))

datatable(best_model, filter = 'top', 
          options = list(pageLength = 7))

Terlihat bahwa model paling terbaik adalah model tentatif tanpa kombinasi drift. Sangat disayangkan, dari 4 kombinasi drift yang ada, hanya kombinasi 1 yang memiliki kategori model terbaik. Sehingga model yang akan digunakan adalah model_nvda.da.

4.4 Overfitting

Overfitting dilakukan dilakukan dengan menaikkan orde AR(p) dan MA(q) dari model terpilih ARIMA(0,1,1) untuk melihat apakah terdapat model lain yang lebih baik dari model saat ini. Kandidat model overfitting adalah ARIMA(1,1,1) dan ARIMA(0,1,2).

Perlu diingat bahwa sebelumnya saya menjelaskan bahwa fungsi model_tentatif yang kami buat itu sudah sekaligus overfitting, dan itu juga sudah membandingkan mana model terbaiknya. Sehingga sebenarnya tidak perlu dilakukan overfitting. Namun saya akan tunjukkan kembali.

datatable(model.tentaif_nvda, filter = 'top', 
          options = list(pageLength = 2))

Terlihat bahwa hanya terdapat model ARIMA(0,1,2) saja. Itu dikarenakan model ARIMA(1,1,1) tidak tersedia di Tabel EACF, sehingga perlu dilakukan manual.

Arima(train_nvda.diff, order=c(1,1,1), method="ML") %>% summary() 
## Series: train_nvda.diff 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.0040  -0.9934
## s.e.  0.0357   0.0081
## 
## sigma^2 = 9.879e-07:  log likelihood = 4460.68
## AIC=-8915.37   AICc=-8915.34   BIC=-8901.27
## 
## Training set error measures:
##                         ME         RMSE          MAE MPE MAPE      MASE
## Training set -6.096225e-06 0.0009921142 0.0006814358 Inf  Inf 0.7812342
##                      ACF1
## Training set -0.001297599

memiliki AIC yang lebih besar dan tidak semua parameternya signifikan. Sehingga tetap model terbaiknya adalah ARIMA(0,1,1).

arima.best_model_nvda <- model_nvda.da

arma.order <- c(
  arima.best_model_nvda[["call"]][["order"]][[2]],
  arima.best_model_nvda[["call"]][["order"]][[4]]
)

5 Analisis Sisaan

Model terbaik hasil identifikasi kemudian dicek asumsi sisaannya. Sisaan model ARIMA harus memenuhi asumsi normalitas, kebebasan sisaan, dan kehomogenan ragam. Diagnostik model dilakukan secara eksplorasi dan uji formal.

5.1 Eksplorasi Sisaan

#Eksplorasi 
sisaan_nvda.da <- model_nvda.da$residuals 
par(mfrow=c(2,2)) 
# QQ-Plot
qqnorm(sisaan_nvda.da) 
qqline(sisaan_nvda.da, col = "dodgerblue3", lwd = 2.5) 
# Plot sisaan
plot(c(1:length(sisaan_nvda.da)), sisaan_nvda.da, 
     main = "Plot Sisaan Model ARIMA",
     xlab = "Periode",
     ylab = "Nilai Sisaan")
abline(h = 0, col = "firebrick2", lty = 2, lwd=2.5)
#Plot ACF dan PACF
stats::acf(sisaan_nvda.da, ylim=c(-.1,.1)) 
stats::pacf(sisaan_nvda.da, ylim=c(-.1,.1)) 

Berdasarkan plot kuantil-kuantil normal, secara eksplorasi ditunjukkan sisaan menyebar tidak normal ditandai dengan tititk-titiknya cenderung tidak mengikuti garis \(45^{\circ}\). Kemudian dapat dilihat juga lebar pita sisaan yang cenderung sama menandakan bahwa sisaan memiliki ragam yang homogen. Plot ACF dan PACF sisaan ARIMA(0,1,1) juga tidak signifikan pada 20 lag awal yang menandakan saling bebas. Kondisi ini akan diuji lebih lanjut dengan uji formal.

5.2 Uji Formal

formal <- data.frame(
  Menguji = c("Kenormalan Siaan", "Kebebasan Sisaan", 
              "Homogenitas Sisaan", "Nilai Tengah = 0"),
  Uji = c("Kolmogorov-Smirnov", "Ljung-Box", 
          "Ljung-Box Sisaan^2", "Uji-T"),
  p.value = c(ks.test(sisaan_nvda.da, "pnorm")[["p.value"]],
              Box.test(sisaan_nvda.da, type = "Ljung")[["p.value"]],
              Box.test((sisaan_nvda.da)^2, type = "Ljung")[["p.value"]],
              t.test(sisaan_nvda.da, mu=0, conf.level=.95)[["p.value"]]) 
)
formal$Keputusan[1] <- ifelse(formal$p.value[1] < .05, "Sisaan Tidak Menyebar Normal",
                              "Sisaan Menyebar Normal")
formal$Keputusan[2] <- ifelse(formal$p.value[2] < .05, "Sisaan Tidak Saling Bebas",
                              "Sisaan Saling Bebas")
formal$Keputusan[3] <- ifelse(formal$p.value[3] < .05, "Ragam Sisaan Tidak Homogen",
                              "Ragam Sisaan Homogen")
formal$Keputusan[4] <- ifelse(formal$p.value[4] < .05, "Nilai Tengah Sisaan \u2260 0",
                              "Nilai Tengah Sisaan = 0")

formal$p.value <- ifelse(round(formal$p.value, 2) == 0, 
                                  "< 2.2e-16",
                                  round(formal$p.value, 2) )

datatable(formal, options = list(pageLength = 4))

5.2.1 Sisaan Menyebar Normal

Selain dengan eksplorasi, asumsi tersebut dapat diuji menggunakan uji formal. Pada tahapan ini uji formal yang digunakan untuk normalitas adalah uji Kolmogorov-Smirnov (KS). Hipotesis pada uji KS adalah sebagai berikut.

\(H_0\) : Sisaan menyebar normal

\(H_1\) : Sisaan tidak menyebar normal

ks.test(sisaan_nvda.da, "pnorm")
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan_nvda.da
## D = 0.49834, p-value < 2.2e-16
## alternative hypothesis: two-sided
if (ks.test(sisaan_nvda.da, "pnorm")[["p.value"]] < 0.05) {
  cat("Karena p-value < 0.05, Maka Sisaan Tidak Menyebar Normal")
} else {
  cat("Karena p-value > 0.05, Maka Sisaan Menyebar Normal")
}
## Karena p-value < 0.05, Maka Sisaan Tidak Menyebar Normal

Berdasarkan uji KS tersebut, didapat p-value yang sangat kecil yakni kurang dari \(2.2\times10^{-16}\). p-value nya kurang dari taraf nyata \(5\%\) sehingga tolak \(H_0\) dan menandakan bahwa sisaan tidak menyebar normal. Sesuai dengan qqplot sisaan.

5.2.2 Sisaan saling bebas

Selanjutnya akan dilakukan uji formal untuk kebebasan sisaan menggunakan uji Ljung-Box. Hipotesis yang digunakan adalah sebagai berikut.

\(H_0\) : Sisaan saling bebas

\(H_1\) : Sisaan tidak tidak saling bebas

Box.test(sisaan_nvda.da, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  sisaan_nvda.da
## X-squared = 0.005934, df = 1, p-value = 0.9386
if (Box.test(sisaan_nvda.da, type = "Ljung")[["p.value"]] < 0.05) {
  cat("Karena p-value < 0.05, Maka Sisaan Tidak Saling Bebas")
} else {
  cat("Karena p-value > 0.05, Maka Sisaan Saling Bebas")
}
## Karena p-value > 0.05, Maka Sisaan Saling Bebas

Berdasarkan uji Ljung-Box tersebut, p-value nya lebih besar dari taraf nyata \(5\%\) sehingga tak tolak \(H_0\) dan menandakan bahwa sisaan saling bebas. Yang berarti tidak ada autokorelasi. Hal ini sesuai dengan hasil eksplorasi menggunakan plot ACF dan PACF.

5.2.3 Sisaan homogen

Hipotesis yang digunakan untuk uji kehomogenan ragam adalah sebagai berikut.

\(H_0\) : Ragam sisaan homogen

\(H_1\) : Ragam sisaan tidak homogen

Box.test((sisaan_nvda.da)^2, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  (sisaan_nvda.da)^2
## X-squared = 4.2891, df = 1, p-value = 0.03836
if (Box.test((sisaan_nvda.da)^2, type = "Ljung")[["p.value"]] < 0.05) {
  cat("Karena p-value < 0.05, Maka Sisaan Tidak Homogen")
} else {
  cat("Karena p-value > 0.05, Maka Sisaan Homogen")
}
## Karena p-value < 0.05, Maka Sisaan Tidak Homogen

Berdasarkan uji Ljung-Box terhadap sisaan kuadrat tersebut, p-value nya lebih kecil dari taraf nyata \(5\%\) sehingga tolak \(H_0\) dan menandakan bahwa ragam sisaan tidak homogen.

5.2.4 Nilai tengah sama dengan nol

Terakhir, dengan uji-t, akan dicek apakah nilai tengah sisaan sama dengan nol. Hipotesis yang diujikan sebagai berikut.

\(H_0\) : nilai tengah sisaan sama dengan 0

\(H_1\) : nilai tengah sisaan tidak sama dengan 0

t.test(sisaan_nvda.da, mu=0, conf.level=.95)
## 
##  One Sample t-test
## 
## data:  sisaan_nvda.da
## t = -0.18269, df = 812, p-value = 0.8551
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -7.470093e-05  6.198009e-05
## sample estimates:
##     mean of x 
## -6.360421e-06
if (t.test(sisaan_nvda.da, mu=0, conf.level=.95)[["p.value"]] < 0.05) {
  cat("Karena p-value < 0.05, Maka Nilai Tengah \u2260 0")
} else {
  cat("Karena p-value > 0.05, Maka Nilai Tengah = 0")
}
## Karena p-value > 0.05, Maka Nilai Tengah = 0

Berdasarkan uji-t tersebut, p-value lebih besar dari taraf nyata \(5\%\) sehingga tak tolak \(H_0\) dan menandakan bahwa nilai tengah sisaan sama dengan nol.

6 Model ARCH-GARCH

Referensi : Measuring Volatility Using the ARCH Model, Volatility Modeling with R : ARCH and GARCH Models

6.1 ARCH

ARCH merupakan singkatan dari Auto-Regressive Conditional Teteroskedasticity.

  1. AutoRegressive (AR): Ini merujuk pada ide bahwa variansi (volatilitas) suatu variabel bergantung pada nilai-nilai sebelumnya dalam deret waktu. Dengan kata lain, variabilitas saat ini dipengaruhi oleh nilai-nilai sebelumnya.
  2. Conditional (Conditional): Ini menunjukkan bahwa perhitungan volatilitas dilakukan secara kondisional terhadap informasi sebelumnya dalam deret waktu. Dengan kata lain, variabilitas diukur dengan mempertimbangkan informasi kondisional terkini.
  3. Heteroskedasticity: Ini merujuk pada sifat ketidakseragaman (Heterogen/Tidak Homogen) dalam variabilitas. Dalam konteks ARCH, heteroskedasticity menunjukkan bahwa variabilitas suatu variabel tidak konstan, tetapi bervariasi sepanjang waktu berdasarkan informasi kondisional.

Jadi, AutoRegressive Conditional Heteroskedasticity (ARCH) adalah model statistik yang digunakan untuk mengukur dan memodelkan volatilitas dalam deret waktu finansial atau ekonometrik. Model ini mengasumsikan bahwa variabilitas harga atau hasil suatu variabel dapat bervariasi seiring waktu dan dipengaruhi oleh nilai-nilai sebelumnya dalam deret waktu.

Volatilitas sendiri merujuk pada tingkat fluktuasi atau variasi harga atau hasil suatu variabel dalam periode waktu tertentu. Volatilitas tinggi menunjukkan fluktuasi yang besar dalam nilai-nilai variabel tersebut, sementara volatilitas rendah menunjukkan fluktuasi yang lebih terbatas. Dalam konteks finansial, volatilitas sering kali dianggap sebagai ukuran risiko, karena fluktuasi harga yang besar dapat menunjukkan ketidakpastian atau risiko yang tinggi.

Model ARCH adalah model yang populer karena spesifikasi variansnya dapat menangkap fitur umum dari deret waktu variabel keuangan. Model ini berguna untuk memodelkan volatilitas dan khususnya perubahan volatilitas seiring waktu. Mari kita lihat kembali pada plot NVDA untuk mengevaluasi volatilitasnya.

#Plot
plot.ts(train_nvda.diff, lty=1, xlab="Periode (Hari)", 
        col = col.nvda[1], lwd = 2,
        main="NVDA", cex.main = 2)

Seperti yang dapat Anda lihat, nilai dari deret ini berubah dengan cepat dari periode ke periode dalam suatu cara yang tampaknya tidak dapat diprediksi, yang dapat kita sebut sebagai volatil. Selain itu, ada periode-periode ketika perubahan besar diikuti oleh perubahan besar berikutnya dan periode di mana perubahan kecil diikuti oleh perubahan kecil berikutnya, merupakan bukti lain dari volatilitas.

6.1.1 Uji Efek ARCH

Uji Langrange multiplier (LM) sering digunakan untuk menguji keberadaan efek ARCH.

Untuk melakukan uji efek ARCH, kita harus

  1. Mengestimasi persamaan rata-rata yang, dalam contoh ini, adalah \(r_t=\beta_0 + e_t\), dimana \(r\) = train_nvda.diff

  2. Mendapatkan residual yang diestimasi \(e^t\)

  3. Mengestimasi \(\hat{e}^2_t = \gamma_0 + \gamma_1 \hat{e}^2_{t-1} + vt\)
    dimana \(vt\) adalah suatu error acak.

Hipotesis

\(H_0\) : \(\gamma_1=0\) (Tidak ada efek ARCH)

\(H_1\) : \(\gamma_1 \ne0\) (Ada efek ARCH)

install_load('dynlm')
# Langkah 1: Mengestimasi persamaan rata-rata r = beta + error
byd.mean <- dynlm(train_nvda.diff ~ 1)

# Langkah 2: Mengambil residual dari model sebelumnya dan mengkuadratkannya
ehatsq <- ts(resid(byd.mean)^2)

# Langkah 3: Regresikan residual kuadrat pada residual kuadrat yang tertinggal satu lag
byd.arch <- dynlm(ehatsq ~ L(ehatsq), data = ehatsq)

summary(byd.arch)
## 
## Time series regression with "ts" data:
## Start = 2, End = 813
## 
## Call:
## dynlm(formula = ehatsq ~ L(ehatsq), data = ehatsq)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -2.235e-06 -9.005e-07 -7.837e-07 -2.000e-07  2.875e-05 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 9.114e-07  8.638e-08  10.552   <2e-16 ***
## L(ehatsq)   7.268e-02  3.504e-02   2.074   0.0384 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.257e-06 on 810 degrees of freedom
## Multiple R-squared:  0.005282,   Adjusted R-squared:  0.004054 
## F-statistic: 4.301 on 1 and 810 DF,  p-value: 0.0384

Kita dapat memeriksa efek ARCH dengan menggunakan fungsi ArchTest() dari paket FinTS. Kita akan menggunakan tingkat signifikansi \(\alpha=0.05\) untuk uji hipotesis nol kita. Perlu diperhatikan bahwa ArchTest() memiliki default lags=12.

install_load('FinTS')
lag <- 12
arch.effect <- data.frame(p.value = numeric(0), 
                          Keterangan = character(0))
for(i in 1:lag){
  p.value <- ArchTest(sisaan_nvda.da, lags=i, demean=TRUE)[["p.value"]]
  arch.effect <- rbind(arch.effect, 
      data.frame(p.value = ifelse(round(p.value, 2) == 0, 
                                  sprintf("%.2e", p.value),
                                  round(p.value, 2) ),
                 Keterangan = ifelse(p.value < 0.05, 
                                     "Signifikan", "Tidak signifikan")
                 ) )
  rownames(arch.effect)[i] <- paste0("Lag ", i)
}
datatable(arch.effect, filter = 'top', 
          options = list(pageLength = 6))

Dari Lag 1 sampai Lag 12 yang ditest, Terlihat bahwa Lag 2 memiliki p-value sebesar \(0.08\) yang lebih besar dari taraf nyata \(5\%\), sehingga tolak \(H_0\). Artinya Lag 2 Tidak sinifikan atau tidak ada efek ARCH(2). Selain Lag 2, Lag nya signifikan. Karena pada saat Lag 2 saja sudah tidak signifikan maka tidak perlu dilakukan model GARCH.

6.1.2 Model ARCH

install_load('rugarch')

create_spec <- function(i, j, arma_order, mean=mean) {
  return(ugarchspec(
    variance.model = list(model = "sGARCH", garchOrder = c(i, j)),
    mean.model = list(armaOrder = arma_order, 
              #model tidak harus memasukkan komponen rata-rata (intercept)
                      include.mean = mean)
  ))
}


arch.garch <- function(data, pp, qq = 0, arma.order = c(0, 0), mean=FALSE) {
  arch.model <- data.frame(Model = character(0),
                            Signif = numeric(0),
                            AIC = numeric(0),
                            Sig.par = character(0),
                            Non_Sig.par = character(0))
  
  best_model <- NULL
  best_sig <- -Inf
  best_aic <- Inf
  
  if(qq == 0){
    n <- 0
  } else{
    n <- 1
  }
  for (i in 1:pp) {
    for (j in n:qq) {
        # GARCH model
        formula <- as.formula(paste("~ garch(", i, ",", j, ")", sep = ""))
        archSpec <- create_spec(i, j, arma.order, mean=mean)
  
      archFit <- ugarchfit(spec = archSpec, data = data, solver = "hybrid")
  
      if (anyNA(coef(archFit)) || anyNA(archFit@fit[["tval"]])) {
        break
      }
  
      aic <- infocriteria(archFit)[[1]] # AIC
      p_values <- archFit@fit$robust.matcoef[, 4] # P-value
      
      # persentase banyaknya param yang signif
      sig <- sum(p_values < 0.05) / length(p_values)
  
      sig_param <- names(p_values[p_values < 0.05])
      non_sig_param <- names(p_values[p_values >= 0.05])
  
      model <- if (qq == 0) {
        paste0("ARCH(", i, ")")
      } else {
        paste0("GARCH(", i, ",", j, ")")
      }
  
      arch.model <- rbind(
        arch.model,
        data.frame(
          Model = model,
          Signif = sig,
          AIC = aic,
          Sig.par = paste(sig_param, collapse = ", "),
          Non_Sig.par = paste(non_sig_param, collapse = ", ")
        )
      )
  
      # Identifikasi best model
      if (sig > best_sig) {
        if (aic < best_aic) {
          best_model <- model
          best_sig <- sig
          best_aic <- aic
        }
      }
    }
  }
  
  return(list(arch.model = arch.model, best_model = best_model))
}

Model ARCH tanpa Intersep

model.arch1 <- arch.garch(data = train_nvda.diff, pp = 12, qq = 0, arma.order = arma.order)
datatable(model.arch1$arch.model, filter = 'top', 
          options = list(pageLength = 6))
cat("Best Model:", model.arch1$best_model)
## Best Model: ARCH(1)

Setelah mencoba iterasi lag 1 sampai lag 12, terlihat bahwa hanya ada satu model yang semua parameternya signifikan yakni ARCH(1). Karena itu tidak perlu membandingkan nilai AIC. Sehingga model terbaik adalah ARIMA(0,1,1)+ARCH(1).

Model ARCH dengan intersep

model.arch2 <- arch.garch(data=train_nvda.diff, pp=12, qq=0, arma.order=arma.order,
                         mean=TRUE)
datatable(model.arch2$arch.model, filter = 'top', 
          options = list(pageLength = 6))
cat("Best Model:", model.arch2$best_model)
## Best Model: ARCH(1)

Sama seperti model ARCH tanpa intersep, model terbaiknya adalah ARCH(1).

#Model dengan intersep
archSpec <- ugarchspec(
  variance.model=list(model="sGARCH", garchOrder=c(1,0)),
  mean.model=list(armaOrder = arma.order, include.mean = TRUE) )
model.arch_nvda1 <- ugarchfit(spec=archSpec, data=train_nvda.diff, solver="hybrid")
coef(model.arch_nvda1)
##            mu           ma1         omega        alpha1 
## -5.087292e-05  4.308273e-02  7.454659e-07  2.989799e-01
#Model tanpa intersep
archSpec <- ugarchspec(
  variance.model=list(model="sGARCH", garchOrder=c(1,0)),
  mean.model=list(armaOrder = arma.order, include.mean = FALSE) )
model.arch_nvda2 <- ugarchfit(spec=archSpec, data=train_nvda.diff, solver="hybrid")
coef(model.arch_nvda2)
##          ma1        omega       alpha1 
## 4.514680e-02 7.495716e-07 2.961780e-01

Model yang dihasilkan adalah seperti diatas.

6.1.3 Plot

plot.ts(ts(model.arch_nvda1@fit$residuals), col=col.nvda[1], lty=1, 
        xlab="Periode (Hari)", ylab="Conditional Variance",
        lwd = 2, main="ARIMA(0,1,1)-ARCH(1) NVDA", cex.main = 1.5)

Terlihat bahwa ada periode-periode yang volatilitas nya tinggi.

6.2 GARCH

Referensi : GARCH models with R programming : a practical example with TESLA stock

Ya walaupun diatas sudah dijelaskan “tidak usah melanjutkan ke garch”, tapi kami hanya ingin mencoba. “Apa yang terjadi jika kita tetap melakukan model GARCH?”. Model terbaik dari GARCH(p,q) ditentukan dari banyaknya paramter yang signifikan dan AIC nya yang terkecil.

model.garch <- arch.garch(data = train_nvda.diff, pp = 5, qq = 5, arma.order = arma.order)
datatable(model.garch$arch.model, filter = 'top', 
          options = list(pageLength = 6))
cat("Best Model:", model.garch$best_model)
## Best Model: GARCH(1,1)

Terlihat bahwa memang data tidak cocok untuk model GARCH(p,q) karena dari rentang lag 1 sampai 12 yang dihasilkan model terbaiknya adalah GARCH(12,12). Namun jika pada modelnya, hanya persentase signifikannya sangat kecil. Sehingga tidak mungin untuk dijadikan model terbaik.

6.3 Model Terbaik

Sehingga diperoleh model terbaik yakni ARIMA(0,1,1)+ARCH(1).

plot(model.arch_nvda1, which="all")
## 
## please wait...calculating quantiles...

Diperoleh mean model yaitu ARIMA(0,1,1) dan varian model yaitu ARCH(1) untuk digunakan dalam peralaman, pada uji dignostik model pada plot QQ-norm terlihat bahwa sebaran residual tidak mengikuti garis lurus sehingga dapat disimpulkan bahwa residual tidak menyebar normal.

Box.test((sisaan_nvda.da)^2, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  (sisaan_nvda.da)^2
## X-squared = 4.2891, df = 1, p-value = 0.03836

7 Akurasi

Akan dilakukan peramalan sebanyak data uji untuk mengukur akurasi. Peramalan dilakukan menggunakan fungsi forecast() . Contoh peramalan berikut ini dilakukan untuk \(h\) hari ke depan. Dengan \(h\) adalah banyaknya data test per perusahaan.

#--Forecast--
  #ARIMA
  ramalan_nvda.da <- forecast::forecast(model_nvda.da, h = nrow(test_nvda)) 
  data.ramalan_nvda.da <- ramalan_nvda.da$mean
  
  #ARCH + Intersep
  ramalan_nvda.arch1 = ugarchforecast(model.arch_nvda1, data = train_nvda.diff, 
                        n.ahead = nrow(test_nvda))
  data.ramal_nvda.arch1 <- ramalan_nvda.arch1@forecast[["seriesFor"]]
  
  #ARCH - Intersep
  ramalan_nvda.arch2 = ugarchforecast(model.arch_nvda2, data = train_nvda.diff, 
                        n.ahead = nrow(test_nvda))
  data.ramal_nvda.arch2 <- ramalan_nvda.arch2@forecast[["seriesFor"]]

7.1 Plot Ramalan

Model ARIMA(0,1,1)

#Plot
plot(ramalan_nvda.da, main="NVDA", xlab="Periode(Hari)", cex.main=4, 
     lwd=2.5, col = col.nvda[1])

Berdasarkan hasil plot ramalan di atas, dapat dilihat bahwa semua ramalan cenderung stabil hingga akhir periode.

Model ARIMA(0,1,1)+ARCH(1)

plot(ramalan_nvda.arch1, which= 1)

Berdasarkan hasil plot ramalan di atas, dapat dilihat bahwa semua ramalan cenderung stabil hingga akhir periode.

#--Diff inv--
  #ARIMA
  ramal_nvda <- diffinv(data.ramalan_nvda.da, differences=1) +
    train_nvda.ts.new[length(train_nvda.ts)] #nilai akhir data latih
  
  #ARCH + Intersep
  ramal_nvda.arch1 <- diffinv(data.ramal_nvda.arch1, differences=1) +
    train_nvda.ts.new[length(train_nvda.ts)] #nilai akhir data latih
  
  #ARCH - Intersep
  ramal_nvda.arch2 <- diffinv(data.ramal_nvda.arch2, differences=1) +
    train_nvda.ts.new[length(train_nvda.ts)] #nilai akhir data latih

#--Transformasi Balik--
  #ARIMA
  ramal_nvda <- 1/(ramal_nvda)^2
  
  #ARCH + Intersep
  ramal_nvda.arch1 <- 1/(ramal_nvda.arch1)^2

  #ARCH - Intersep
  ramal_nvda.arch2 <- 1/(ramal_nvda.arch2)^2

#Data Frame
ramal <- data.frame(
  Date = data2 %>% dplyr::select(Date)  ,
  `ARIMA` = c(as.vector(train_nvda.ts), ramal_nvda[-1] ),
  `ARCH+Intersep` = c(as.vector(train_nvda.ts), ramal_nvda.arch1[-1] ),
  `ARCH-Intersep` = c(as.vector(train_nvda.ts), ramal_nvda.arch2[-1] ),
  Name = rep(c("NVDA"), 
              each = nrow(train_nvda) + nrow(test_nvda) )
)
colnames(ramal) <- c("Date", "ARIMA", "ARCH+Intersep", "ARCH-Intersep", "Name")
chart <- 
ggplot() +
# Akhir data
  geom_vline( #Buat garis batas data
    xintercept = as.numeric(as.Date( max(data2$Date) )), 
             linetype = "dotted", color = "red") +
  #Label Data Asli 
  annotate( "rect", alpha=.1, fill="#4EC2C1",
            xmin=as.Date(min(data2$Date)), 
            xmax=as.Date(data2$Date[nrow(train_nvda)]),
            ymin=max(nvda$`Adj Close`) * 1.25, ymax=Inf ) + 
  
  annotate( "text", color="#4EC2C1",
            x = as.Date(data2$Date[ceiling(nrow(train_nvda)/2)]), 
            y = max(nvda$`Adj Close`) * 1.3, 
    label = "Data Asli", size=10) + 
  annotate( "text", color="#4EC2C1",
            x = as.Date(data2$Date[ceiling(nrow(train_nvda)/2)]), 
            y = max(nvda$`Adj Close`) * 1.2, 
    label = paste0(nrow(train_nvda)," Hari (", 
                   round(nrow(train_nvda)/nrow(nvda)*100, 1), 
                   "%)"), size=10) + 
  
  #Label Data Ramal
  annotate( "rect", alpha=.1, fill="violetred",
            xmin=as.Date(data2$Date[nrow(train_nvda)]), 
            xmax=as.Date(max(data2$Date)),
            ymin=max(nvda$`Adj Close`) * 1.25, ymax=Inf ) + 
  
  annotate( "text", color="violetred",
            x = as.Date(data2$Date[ceiling( .89*length(data2$Date)/6)]) , 
            y = max(nvda$`Adj Close`) * 1.3, 
    label = "Data Ramal", size=10) +
  annotate( "text", color="violetred",
            x = as.Date(data2$Date[ceiling( .89*length(data2$Date)/6)]) , 
            y = max(nvda$`Adj Close`) * 1.2, 
    label = paste0(nrow(test_nvda)," Hari (", 
                   round(nrow(test_nvda)/nrow(nvda)*100, 1), 
                   "%)"), size=10) +
  
  geom_vline( #Buat garis batas data
    xintercept = as.Date(data2$Date[nrow(train_nvda)]) , 
             linetype = "dotted", color = "black", linewidth = 2) +

#Time Series
  #Data RAMAL ARIMA
    #Dari data asli & ramal
  geom_line(data = ramal, linewidth=2,
            aes(x = Date, y = `ARIMA`, col = "ARIMA")) +
  #Data Ramal ARCH + Intersep
    #Dari data ramal saja
  geom_line(data = ramal[nrow(train_nvda):nrow(nvda),], linewidth=2,
            aes(x = Date, y = `ARCH+Intersep`, col = "ARCH+")) +
  #Data Ramal ARCH - Intersep
    #Dari data ramal saja
  geom_line(data = ramal[nrow(train_nvda):nrow(nvda),], linewidth=2,
            aes(x = Date, y = `ARCH-Intersep`, col = "ARCH-")) +
  #Data Test
  geom_line(data = test_nvda, linewidth=2,
            aes(x = Date, y = `Adj Close`, col = "Test")) +
  
  #Penanda
  geom_point(data = tail(train_nvda, 1), alpha = .5, 
             aes(x = Date, y = `Adj Close`), stroke=2,
             size = 15, shape = 21, color = "black", fill="violetred") +
  scale_colour_manual(values = c("ARIMA"=col.nvda[1],"ARCH+"="#539282",
                                "ARCH-"="#93C6B9","Test"=col.nvda[2]
                                 )) +
  theme.ts + #THeme
  labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
       title = "Akurasi Ramalan Saham NVDA",
       subtitle = "Perbandingan Hasil Ramalan dengan data Asli\n") +
  #Axis
    coord_cartesian(clip = "off"
  ) +
    scale_x_date( #Sumbu x
    date_breaks = "1 year",  # Menampilkan label setiap tahun
    date_labels = "%Y",  # Format label tahun
    limits = c(as.Date(min(ramal$Date)), 
               as.Date(max(ramal$Date)) + 80)
  ) +
    scale_y_continuous( #Sumbu y
    labels = scales::dollar_format(prefix = "$") #tambahin dolar
  ) +
  
#Bagi Data Train & Test
  geom_segment(aes(x = as.Date(max(train_nvda$Date)), 
                   xend = as.Date(max(train_nvda$Date)), 
                   y = max(nvda$`Adj Close`)* .58, 
                   yend = max(nvda$`Adj Close`)*.7), 
               arrow = arrow(type = "closed", length = unit(0.1, "inches")), 
               lineend = "round", color = "gray60", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(train_nvda$Date)), 
                      y = max(nvda$`Adj Close`)*.7, 
                      label = max(train_nvda$Date) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "gray60", box.color = "white", parse = TRUE
  ) +

#Data Uji
  geom_segment(aes(x = as.Date(max(data2$Date)), 
                   xend = as.Date(max(data2$Date))+60, 
                   y = nvda$`Adj Close`[nrow(nvda)], 
                   yend = nvda$`Adj Close`[nrow(nvda)]), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "gray60", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(data2$Date))+60, 
                      y = nvda$`Adj Close`[nrow(nvda)], 
                      label = "Data Asli" ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "gray60", box.color = "white", parse = TRUE
  ) +

# Data Ramal ARIMA
  geom_segment(aes(x = as.Date(max(data2$Date)), 
                   xend = as.Date(max(data2$Date)) + 60, 
                   y = max(ramal$`ARIMA`), 
                   yend = max(ramal$`ARIMA`)), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "#D02A49", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(data2$Date)) +60, 
                      y = max(ramal$`ARIMA`), 
                      label = "ARIMA" ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "#D02A49", box.color = "white", parse = TRUE
  ) +
  
# Data Ramal ARCH + Intersep
  geom_segment(aes(x = as.Date(max(data2$Date)), 
                   xend = as.Date(max(data2$Date)) + 70, 
                   y = max(ramal$`ARCH+Intersep`), 
                   yend = max(ramal$`ARCH+Intersep`)), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "#4B75BA", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(data2$Date)) +70, 
                      y = max(ramal$`ARCH+Intersep`), 
                      label = "ARCH dengan <br> Intersep" ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "#4B75BA", box.color = "white", parse = TRUE
  ) +
  
# Data Ramal ARCH - Intersep
  geom_segment(aes(x = as.Date(max(data2$Date)), 
                   xend = as.Date(max(data2$Date)) + 70, 
                   y = ramal$`ARCH-Intersep`[nrow(nvda)], 
                   yend = ramal$`ARCH-Intersep`[nrow(nvda)]), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "#DC4F26", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(data2$Date)) +70, 
                      y = ramal$`ARCH-Intersep`[nrow(nvda)], 
                      label = "ARCH tanpa <br> Intersep" ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "#DC4F26", box.color = "white", parse = TRUE
  ) +

#Tanggal akhir
  geom_segment(aes(x = as.Date(max(data2$Date)), 
                   xend = as.Date(max(data2$Date)), 
                   y = -Inf, 
                   yend = max(nvda$`Adj Close`)*.15), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "gray60", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(data2$Date)), 
                      y = max(nvda$`Adj Close`)*.15, 
                      label = max(data2$Date) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "gray60", box.color = "white", parse = TRUE
  )

chart

#Export Chart
ggsave("06_Akurasi_model.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 27)

Dapat dilihat bahwa rata-rata harga saham Amazon diramalkan akan cenderung sedikit menurun setiap periodenya. Di sisi lain, saham NVIDIA diramalkan akan menaik pesat setiap periodenya. Sedangkan harga saham lainnya cenderung naik sedikit.

7.2 Perbandingan Data Aktual dengan Data Ramal

#Perbandingan
perbandingan_nvda.da <- matrix(data=c(head(test_nvda.ts, n=nrow(test_nvda)), 
                                      ramal_nvda[-1],
                                      ramal_nvda.arch1[-1],
                                      ramal_nvda.arch2[-1]),
                               nrow = nrow(test_nvda), ncol = 4)
colnames(perbandingan_nvda.da) <- c("Aktual","ARIMA", 
                                    "ARCH+Intersep","ARCH-Intersep")
datatable(perbandingan_nvda.da, filter = 'top', 
          options = list(pageLength = 5))

7.3 Akurasi

#Akurasi Model ARIMA
akurasi.nvda <- accuracy(ts(ramal_nvda[-1]), 
                         head(test_nvda.ts, n=nrow(test_nvda) )) %>%
  as.data.frame() %>%
  cbind(Model = ramalan_nvda.da[["method"]]) %>% 
  relocate(Model, .before = 1)
row.names(akurasi.nvda) <- "ARIMA"

#Akurasi Model ARCH + Intersep
akurasi.nvda.arch1 <- accuracy(ts(ramal_nvda.arch1[-1]), 
                         head(test_nvda.ts, n=nrow(test_nvda) )) %>%
  as.data.frame() %>%
  cbind(Model = model.arch1$best_model) %>% 
  relocate(Model, .before = 1)
row.names(akurasi.nvda.arch1) <- "ARCH+Intersep"

#Akurasi Model ARCH - Intersep
akurasi.nvda.arch2 <- accuracy(ts(ramal_nvda.arch2[-1]), 
                         head(test_nvda.ts, n=nrow(test_nvda) )) %>%
  as.data.frame() %>%
  cbind(Model = model.arch2$best_model) %>% 
  relocate(Model, .before = 1)
row.names(akurasi.nvda.arch2) <- "ARCH-Intersep"
akurasi <- rbind(akurasi.nvda, akurasi.nvda.arch1, akurasi.nvda.arch2) %>%
  mutate(Keterangan = case_when(
    MAPE < 10 ~ "Sangat Baik",
    MAPE >= 10 & MAPE <= 20 ~ "Baik",
    MAPE > 20 & MAPE <= 50 ~ "Layak",
    MAPE > 50 ~ "Tidak Akurat"
  )) %>% relocate(Keterangan, Model, MAPE)

datatable(akurasi, filter = 'top', 
          options = list(pageLength = 7))

Dari tabel diatas terlihat bahwa model ARIMA(0,1,1)-ARCH(1) dengan Intersep berkatergori “Baik” dengan nilai MAPE nya sebesar \(15.64\%\) (kurang dari \(20\%\)).

8 Peramalan

8.1 Stasioneritas

8.1.1 Ragam

index <- seq(1:nrow(nvda)) #Beda titik potong
bc =boxcox(nvda.ts~index, lambda=seq(-2, 4, by=.01))
title("NVDA", cex.main = 2)

lambda_nvda.ramal <- bc$x[which.max(bc$y)] #Nilai Rounded Lambda
sk_nvda.ramal <- bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)] #SK
lsk.ramal <- data.frame(
  Lambda = c(lambda_nvda.ramal ),
  Batas.Bawah = c( min(sk_nvda.ramal)  ),
  Batas.Atas = c(max(sk_nvda.ramal) )
) %>%
  mutate(across(where(is.numeric), ~round(., 2)))

lsk.ramal$Keterangan <- apply(lsk.ramal, 1, function(row) {
  ifelse(row["Batas.Bawah"] <= 1 && row["Batas.Atas"] >= 1,
         "Ragam Stasioner",  "Ragam Tidak Stasioner")
})

rownames(lsk.ramal) <- c("NVDA")

datatable(lsk.ramal, filter = "none")

Data tidak Stasioner dalam ragam.

tran.ramal <- lsk.ramal %>%
  mutate(
    Transformasi = trans(Lambda),
    `Transformasi Balik` = trans_balik(Lambda)
  ) %>% dplyr::select(Lambda, Transformasi, `Transformasi Balik`)

datatable(tran.ramal, filter = 'none')

Berikut transformasi yang cocok.

8.1.2 Transformasi

#Transformasi
nvda.ts.new <- 1/sqrt(nvda.ts)

#Plot BoxCox
index <- seq(1:nrow(nvda))
bc <-boxcox(nvda.ts.new ~index, lambda=seq(-2, 4, by=.01))
title("NVDA", cex.main=2)

lambda_nvda.ramal.new <- bc$x[which.max(bc$y)] #Nilai Rounded Lambda
sk_nvda.ramal.new <- bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)] #SK
lsk.ramal.new <- data.frame(
  Lambda = c(lambda_nvda.ramal.new ),
  Batas.Bawah = c( min(sk_nvda.ramal.new)  ),
  Batas.Atas = c( max(sk_nvda.ramal.new)  )
) %>%
  mutate(across(where(is.numeric), ~round(., 2)))

lsk.ramal.new$Keterangan <- apply(lsk.ramal.new, 1, function(row) {
  ifelse(row["Batas.Bawah"] <= 1 && row["Batas.Atas"] >= 1,
         "Ragam Stasioner",  "Ragam Tidak Stasioner")
})

rownames(lsk.ramal.new) <- c("NVDA")

datatable(lsk.ramal.new, filter = 'none')

Terlihat bahwa data sudah stasioner dalam ragam.

8.1.3 Rataan & Differencing

#Diff
nvda.diff <- diff(nvda.ts.new, differences = 1) 

tseries::adf.test(nvda.diff)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  nvda.diff
## Dickey-Fuller = -8.9643, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
if (tseries::adf.test(nvda.diff)[["p.value"]] < 0.05) {
  cat("Karena p-value < 0.05, Maka Rataan Stasioner")
} else {
  cat("Karena p-value > 0.05, Maka Rataan Tidak Stasioner")
}
## Karena p-value < 0.05, Maka Rataan Stasioner

Setelah di differencing, data sudah stasioner.

8.2 Model

#Model
archSpec <- ugarchspec(
  variance.model=list(model="sGARCH", garchOrder=c(1,0)), #ARCH(1)
  mean.model=list(armaOrder = arma.order, include.mean = TRUE) ) #ARIMA(0,1,1) + Intersep
model.ramal <- ugarchfit(spec=archSpec, data=nvda.diff, solver="hybrid")

#--Forecast--
nDate <- seq(from = max(data2$Date) + 1, 
             to = as.Date("2023-12-31"), #Dampe akhir tahun
             by = "days")

ramalan_nvda <- ugarchforecast(model.ramal, data = nvda.diff, 
                      n.ahead = length(nDate))
data.ramal_nvda <- ramalan_nvda@forecast[["seriesFor"]]

data.ramal_nvda.bb <- ramalan_nvda@forecast[["seriesFor"]]/
  -( ramalan_nvda@forecast[["sigmaFor"]] )^2

data.ramal_nvda.ba <- ramalan_nvda@forecast[["seriesFor"]]/
  ( ramalan_nvda@forecast[["sigmaFor"]] )^2
#--Diff inv--
  #Mean
  ramal_nvda <- diffinv(data.ramal_nvda, differences=1) +
    nvda.ts.new[length(nvda.ts)] #nilai akhir data latih
  #BB
  ramal_nvda.bb <- diffinv(data.ramal_nvda.bb, differences=1) +
    nvda.ts.new[length(nvda.ts)] #nilai akhir data latih
  #BA
  ramal_nvda.ba <- diffinv(data.ramal_nvda.ba, differences=1) +
    nvda.ts.new[length(nvda.ts)] #nilai akhir data latih


#--Transformasi Balik--
  ramal_nvda <- 1/(ramal_nvda)^2
  ramal_nvda.bb <- 1/(ramal_nvda.bb)^2
  ramal_nvda.ba <- 1/(ramal_nvda.ba)^2

#Data Frame
data_ramal <- data.frame(
  Date = seq(from = min(data2$Date),
             to = max(nDate), by="days")  ,
  `Adj Close` = c(as.vector(nvda.ts), ramal_nvda[-1] ),
  `BB` = c(as.vector(nvda.ts), ramal_nvda.bb[-1] ),
  `BA` = c(as.vector(nvda.ts), ramal_nvda.ba[-1] ),
  Name = rep(c("NVDA"), 
              each = length(nvda.ts) + length(nDate) )
)
colnames(data_ramal) <- c("Date", "Adj Close", "BB", "BA", "Name")

8.3 Plot

harga_akhir <- nvda$`Adj Close`[nrow(nvda)]
harga_ramal <- data_ramal$`Adj Close`[nrow(data_ramal)]
perc <- (harga_ramal-harga_akhir)/harga_akhir * 100

chart <- 
ggplot() +
# Akhir data
  geom_vline( #Buat garis batas data
    xintercept = as.numeric(as.Date( max(data_ramal$Date) )), 
             linetype = "dotted", color = "red") +
  #Label Data Asli 
  annotate( "rect", alpha=.1, fill="#4EC2C1",
            xmin=as.Date("2022-06-01"), 
            xmax=as.Date(nvda$Date[nrow(nvda)]),
            ymin=max(nvda$`Adj Close`) * 1.25, ymax=Inf ) + 
  
  annotate( "text", color="#4EC2C1",
            x = as.Date(data_ramal$Date[ceiling(nrow(data_ramal)*.7)]), 
            y = max(nvda$`Adj Close`) * 1.3, 
    label = "Saham Historis", size=10) + 
  annotate( "text", color="#4EC2C1",
            x = as.Date(data_ramal$Date[ceiling(nrow(data_ramal)*.7)]), 
            y = max(nvda$`Adj Close`) * 1.2, 
    label = paste0(" "), size=10) + 
  
  #Label Data Ramal
  annotate( "rect", alpha=.1, fill="violetred",
            xmin=as.Date(data_ramal$Date[nrow(nvda)]), 
            xmax=as.Date(max(data_ramal$Date)),
            ymin=-Inf, ymax=Inf ) + 
  annotate( "text", color="violetred",
            x = as.Date(max(nDate)-25) , 
            y = min(data2$`Adj Close`) * 10, 
    label = "Data\nRamal", size=10) +
  
  geom_vline( #Buat garis batas data
    xintercept = as.Date(data_ramal$Date[nrow(nvda)]) , 
             linetype = "dotted", color = "black", linewidth = 2) +

#Time Series
  #Data RAMAL 
  geom_line(data = data_ramal, linewidth=2,
            aes(x = Date, y = `Adj Close`, col = "NVDA")) +
  
  #Penanda
  geom_point(data = tail(nvda, 1), alpha = .5, 
             aes(x = Date, y = `Adj Close`), stroke=2,
             size = 15, shape = 21, color = "black", fill="violetred") +
  scale_colour_manual(values = cols) +
  theme.ts + #THeme
  labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
       title = "Ramalan Saham NVDA",
       subtitle = "Peramalan Saham NVDA Hingga Akhir Tahun 2023\n") +
  #Axis
    coord_cartesian(clip = "off"
  ) +
    scale_x_date( #Sumbu x
    date_breaks = "1 year",  # Menampilkan label setiap tahun
    date_labels = "%Y",  # Format label tahun
    limits = c(as.Date("2022-06-01"), 
               as.Date(max(data_ramal$Date)) + 40)
  ) +
    scale_y_continuous( #Sumbu y
    labels = scales::dollar_format(prefix = "$") #tambahin dolar
  ) +
  
#Tanggal akhir data asli
  geom_richtext(
    data = data.frame(x = as.Date(max(nvda$Date)), 
                      y = max(nvda$`Adj Close`)*.75, 
                      label = max(nvda$Date) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "gray60", box.color = "white", parse = TRUE
  ) +
  
# Harga akhir Data Asli
  geom_segment(aes(x = as.Date(max(nvda$Date)), 
                   xend = as.Date(max(nvda$Date))- 30, 
                   y = nvda$`Adj Close`[nrow(nvda)], 
                   yend = max(data_ramal$`Adj Close`)*.9 ) , 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "#4B75BA", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(nvda$Date)) - 30, 
                      y = max(data_ramal$`Adj Close`) *.9, 
                      label = paste0("$ ", harga_akhir %>% round(.,2)) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "#4B75BA", box.color = "white", parse = TRUE
  ) +
  
# Harga akhir Data Ramal
  geom_segment(aes(x = as.Date(max(data_ramal$Date)), 
                   xend = as.Date(max(data_ramal$Date)) + 30, 
                   y = max(data_ramal$`Adj Close`), 
                   yend = max(data_ramal$`Adj Close`)), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "#4B75BA", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(data_ramal$Date)) +30, 
                      y = max(data_ramal$`Adj Close`), 
                      label = paste0("$ ", harga_ramal %>% round(.,2)) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "#4B75BA", box.color = "white", parse = TRUE
  ) +
# Harga akhir data ramal beda brp persen
  annotate( "text", color=ifelse(perc>0, "green4", "red3"),
            x = as.Date(max(data_ramal$Date)) +30, 
            y = max(data_ramal$`Adj Close`) *.94, 
    label = paste0(ifelse(perc>0, "+ ", "- "), perc %>% round(.,2), "%"), 
    size=8 ) + 

#Tanggal akhir ramal
  geom_segment(aes(x = as.Date(max(data_ramal$Date)), 
                   xend = as.Date(max(data_ramal$Date)), 
                   y = -Inf, 
                   yend = max(nvda$`Adj Close`)*.15), 
               arrow = arrow(type = "closed", length = unit(.1, "inches")), 
               lineend = "round", color = "gray60", size=1.5) +
  geom_richtext(
    data = data.frame(x = as.Date(max(data_ramal$Date)), 
                      y = max(nvda$`Adj Close`)*.15, 
                      label = max(data_ramal$Date) ),
    aes(x, y, label = label), size = 8, color = "white", 
    fill = "gray60", box.color = "white", parse = TRUE
  )

chart

#Export Chart
ggsave("07_Ramalan.png", chart, path = export.chart,
        dpi = 300, height = 12, width = 27)

Model ARIMA(0,1,1)-ARCH(1) dengan intersep digunakan untuk meramalkan hingga periode akhir tahun 2023. Pada Gambar 12 terlihat bahwa adjusted close saham NVIDIA diramalkan akan menyentuh angka seharga 555,95 USD atau naik sebesar \(15,02\%\) hingga akhir tahun 2023.

9 Penutup

#Export Data
install_load('openxlsx')
#Model Tentatif 
write.xlsx(list("Tentatif" = model.tentaif_nvda,
                "Drift1" = drift1_nvda,
                "Drift2" = drift2_nvda,
                "Drift3" = drift3_nvda,
                "Drift4" = drift4_nvda,
                "Arch-Intersep" = model.arch1$arch.model, 
                "Arch+Intersep" = model.arch2$arch.model, 
                "Garch" = model.garch$arch.model
                ), 
           file = "Model_NVDA.xlsx")

9.1 Kesimpulan

Kesimpulan yang dapat ditarik pada penelitian ini adalah model ARIMA(0,1,1)-ARCH(1) dengan intersep merupakan model terbaik yang dapat meramalkan adjusted close saham NVIDIA dalam kurung waktu 11 November 2023 hingga 31 Desember 2023. Nilai MAPE yang didapatkan dari model terbaik adalah sebesar \(15.64\%\). Adjusted close saham NVIDIA diramalkan akan naik hingga menyentuh angka \(555.96\) USD atau naik sebesar \(15.02\%\). Hal ini menunjukkan bahwa jika ingin membeli saham harian NVIDIA pada kurung waktu 11 November 2023 hingga 31 Desember 2023 merupakan hal yang tepat dikarenakan saham akan naik.

9.2 Referensi